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FOREWORD 


This report, prepared by the Dynamics and Loads Section, Martin 
Marietta Corpora tion* Denver Division, under Contract NAS8-29946, pre- 
sents the technical approach and the results of a study contract for 
the vibration characteristics of a liquid in a container of arbitrary 
axisymmetric shape with surface tension forces of the same order as 
acceleration forces (Bond Number 1) . The study was administered by 
the National Aeronautics and Space Administration, George C. Marshall 
Space Flight Center, Huntsville, Alabama, under the direction of Mr. 
Frank Bugg, Systems Dynamics Laboratory. 
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1 . INTRODUCTION 


Knowledge of the dynamics of liquid propellant in a low Bond number 
(defined to be mass density times acceleration of gravity and the square 
of a characteristic length divided by the surface tension) environment is 
critical to the design of certain spacecraft systems with respect to or- 
bital propellant transfer and attitude control system. The proposed re- 
usable Space Tug will be required to perform orbital and docking maneu- 
vers with as much as 90% of the vehicle mass as liquid fuel. The effects 
of the liquid mass on the control system will be significant. In the ab- 
sence of normal gravity, liquid free surfaces tend to distort drastically 
under the influence of relatively small disturbances. The propellant mo- 
tions in low gravity may cause inefficient spacecraft operation or even 
mission failure, if the motion has not been accurately accounted for in 
the design of the vehicle. 

Because the Bond number will be near or below unity for a large por- 
tion of the Space Tug flight environment, the equilibrium configuration 
of the liquid and natural frequencies and mode shapes will be primarily 
dependent on liquid surface tension and contact angle. A large body of 
experimental data have been produced on low Bond number sloshing using 
drop tower techniques and at 1-g using very small models,, but there has 
been little successful analytical work applicable for tanks of general 
shape at Bond numbers of one. This study applied the finite element com- 
puter technique to the low Bond number slosh problem for tanks of general 
axisymmetric shape. 

The study resulted in the development of digital computer programs 
for the determination of liquid free surface equilibrium shape, and lat- 
eral slosh natural vibration mode shapes and frequencies for a liquid in 
a container of arbitrary axisymmetric shape with surface tension forces 
the same order of magnitude as acceleration forces (Bond number~l). 

For the vibration analysis, a finite volume element representation of the 
liquid was used. 

The liquid free surface equilibrium shapes were computed for several 
tanks at various contact angles and ullage volumes. One of these config- 
urations was selected for vibration analysis and lateral slosh mode shapes 
and natural frequencies were obtained. 

This report provides documentation of the above results. 



(This page blank) 
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2. TECHNICAL APPROACH 

The problem was approached in two distinct steps: 

a) establish the static equilibrium shape for a given axisymmetric 
container subject to contact angle and ullage volume constraints; 

b) given the static equilibrium shape, model the system as an assem- 
blage of finite elements and obtain natural frequencies and vi- 
bration modes. 

2 , 1 E quilibrium Configuration for a Liquid Bounded by Tree Surface and 
Axisymmetric Container 

Of paramount importance to the solution of the low Bond number slosh 
problem is definition of the equilibrium free surface configuration. We 
seek a solution, for a given container geometry, such that two conditions 
are satisfied. The conditions are: 

a) ullage (or liquid) volume consistent with user specifications, 
and 

b) liquid/container contact angle consistent with user specifica- 
tions . 

In the following sections we develop a formulation for the equilibrium 
surface shape consistent with the above two specifications. The analy- 
tical developments presented herein assume a geometrically axisymmetric 
container subjected to an axisymmetric acceleration field. 

2.1.1 Axisymmetric Meniscus Sha pe - Under the assumptions listed 
above, it is evident that the meniscus configuration will also be axisym- 
metric. Definition of the meniscus shape can be accomplished through 
examination of a force balance. Consider the annular ring cut from the 
meniscus as shown in Figure 2-1. 








s ~ arc length 


5 


and it follows that 

(r+dr)(sin/? cos d/? + cos^ sin d/? ) - r sin /J 

P-P„\ •/ o d£ . n . . dft , - / 0 n 

1 2 1 r (cos /? cos — - sm/? sin—) ds . (2) 

For d/? small we have 
sin d/? — d/? 
cos d/3—1 

and Equation (2) becomes 

(r+dr) (sin/? + cos/? d/? ) - r sin/? 

P 1~ P 2 ^ r (cos/? - sin/? ~ ) ds (3) 


or 

r cos/8 


M + sin o dr „ 
ds P ds 


P -P 
1 2 


r cos/3 


(4) 




where second order terms have been neglected 

Introducing 

. 0 _ dh _ . 1 

s ^p ~— B ~ h 


o dr ~ 1 

cos/S - r 

„ &P d 2 h " 

COS<i Ts = 77 “ h 

ds 


yields 


P.-P, 


or 


r /?* cos/8 + r 1 sin/? - 1 2 ] r cos/? 

O' 


rh ,,; -b r 1 h* = [ III*! \ r r 1 


(5) 

( 6 ) 
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and finally, we have the differential equation for the axisymmetric 
meniscus'* as 


ts (rh,) = 


p -p 
1 2 


rr 


(7) 


A second differential equation arising from the geometry is 
(r 1 ) 2 + (h 1 ) 2 “ 1 ( 8 ) 


a) Li quid Below Meniscus 

For a free meniscus with liquid below (Figure 2-2) we have 

P - P = ullage gas pressure 

1 G 

P = P - P^ - pgh = liquid pressure 

2 L o 


where p * liquid density 
g = acceleration 


P^ - liquid pressure at origin 



* This is exactly the result presented in NASA SP-106, H. N. Abramson, 
editor 
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and the governing differential equations are 
rh 1 + r f h f = rr 1 ^ + 

(r 1 ) 2 + (h') 2 = 1 



b) Liquid Above Meniscus 

For a free meniscus with liquid above (Figure 2-3) we have 
Pf - - Pl o - pgh - liquid pressure 

P T = P = ullage gas pressure 

Li Ct 



Figure 2-3. Meniscus With Liquid Above 


and the governing differential equations are 


rh" + r h = -rr* / G 


i / P r _p L + Pgh 


(r 1 ) 2 + (h') 2 = 1 


( 10 ) 


2.1.2 Volume Contained Within Meniscus - With reference to Figure 
2-1, it is seen that the volume of liquid (or ullage gas) contained within 
the boundary of the meniscus can be represented as 

dV = rr r^dh 


(ID 
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or V' = tt r 2 h * (11) 

and when h‘> 0 then V l > 0 
and when h*< 0 then V*< 0 
so that for 

a. liquid below meniscus V*> 0 is ullage volume 

V'< 0 is liquid volume 

b. liquid above meniscus V*> 0 is liquid volume 

V f < 0 is ullage volume 

2.1.3 The State Equations - Definition of the applicable set of 
state equations follows from Equations (7), (8) and (11). We have 

rh” + rV = ^ ? r*V j rr* 

r'r" + h'h" = 0 (12) 

V' - 7rr 2 h’ 

with P^ = 

P^ = Fl q - pgh for liquid below meniscus, 
and P x =P L C) “ PS h 

P = P for liquid above meniscus. 

2 G 

We now define a set of state variables as 



( 13 ) 


and the governing state equations become 



9 


rh" + r'h' = (A + Bh) rr* 

r'r" + h'h" = 0 (14) 

Wtt = r 2 h* 

with A » )/ ff 

b o 

B = pg/o for liquid below meniscus, 

and A = (?L 

B = -pg/a for liquid above meniscus. 

The system governing state equations follow as 



*2 “ h ' * Y 4 

¥ 3 “ r " = - Y 4 (A + BY 2 ' VV 

Y 4 = h " = Y 3 (A + BY 2 * W . 

Y 5 = Y 1 Y 4 

subject to the initial values 
V x (0) = 0 Y 4 (0) = 0 

Y_ (0) = 0 Y-(0) =0. (16) 

2 6 

y 3 (o) = 1 

The equations can be integrated using a numerical procedure to yield a 
trajectory which defines (for specified values of the coefficients A and 
B) continuous values of h, r and the contained volume V, The apparent 
singularity at the origin occur ing in the 3rd and 4th of Equation (15) 
can be resolved through application of L , Hospital , s Rule. 

2.1.4 Definition of the Container - The development to this point 
has considered only the axisymmetric meniscus shape without regard to the 
nature of the container. For purposes of this investigation, it will be 
convenient to describe the container geometry as shown in Figure 2-4. 
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Figure 2-4. Container Geometry 


2,1.5 Definition of the Contact Angle - Consider the fluid free 
surface intersection with the container. Using the usual definition of 
contact angle (i.e., the angle between the container and the fluid sur- 
face measured through the fluid) we can establish a measure of the con- 
tact angle through examination of Figure 2-5. We define* 


tan y 



x=x a 


i 

r 

c 


tan a 


r'/x' 


x=x a 


dX 


where 


X ~ — 
f ds 


~ = £ (X i + h) = h ' 


* The subscripts f and c refer to fluid and container, respectively; 
the subscript a denotes the fluid/container intersection point. 
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Figure 2-5. Contact Angle Geometry 

and, with reference to the figure, it is evident that the contact angle 
is 

0 = a - y (17) 


2.1.6 Implements tion - The digital computer program to define the 
static equilibrium shape runs in one of two distinct and independent 
modes of operation. The mode selection is controlled by the Boolean 
input control variable SEARCH. If SEARCH is input as .TRUE., the search 
mode is selected; if SEARCH is input as .FALSE., the survey mode is se- 
lected. The same program is used in both cases, however, the input re- 
quirements and the outputs obtained are different enough to justify 
separate discussion of each mode. 

2. 1.6.1 Survey Mode - While developing and checking out 
what is now the search mode, it was found to be quite difficult, in 
many cases, to supply the program with initial values of the iteration 
parameters, A Q and Xj_, sufficiently close to the desired values to guar- 
antee convergence of the solution algorithm. (A complete discussion of 
A c and will be found in Section 2. 1.6. 2.) This difficulty is espe- 
cially severe prior to the user obtaining any computational experience 
with a given container. In other words, the user thinks of the problem 
in terms of contact angle (0) and ullage volume (V u ) but the solution 
algorithm developed from the differential equations derived in Section 
2.1.3 requires A 0 and X± to effect a solution. Therefore, the survey 
mode was developed and incorporated into the program to assist the user 
in generating a mapping, or a transformation from 9, V u coordinates to 
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A , coordinates. It is strongly recommended that the user take ad- 
vantage of this feature before attempting to arrive at a specific solu- 
tion via the search mode. 

The survey mode operates by cycling through a user specified range 
of the two iteration parameters producing a table of 0 and V u vs for 
each A 0 . On option, the table may be plotted, three plot frames (0 vs 
X t> V u vs X is and 0 vs V u ) for each A 0 . This procedure may need to be 
repeated several times with different ranges on A 0 and X^ to obtain the 
desired degree of correspondence but is inexpensive compared to running 
the search mode blind with trial and error initial values. Given this 
survey data and a desired 0 and V u , the user can interpolate and/or cross- 
plot to obtain good initial values of the iteration parameters for input 
to the search mode. 

Specific input requirements for the survey mode are given in Section 
6. 2. 1.1. A sample of the tabular printed output is shown in Figure 3-1, 
and samples of the optional plotted output are displayed in Figures 3-2 
through' 3- 5 . 


2. 1.6.2 Search Mode - The solution algorithm as encoded in 
the search mode of the digital computer program consists of three phases: 
input and initialization, iteration, and output of solution. 

The input and initialization phase reads input data in NAMELIST 
format (as fully described in Section 6.2.1), performs several checks 
on the consistency of the input data, computes required constants, and 
initializes certain variables. 

The iteration phase performs the computations required to arrive 
at the coordinates of the free surface static equilibrium shape. The 
four nested loops which comprise the iteration phase logic are, from 
inner to outer, the integration loop, the X loop, the A loop, and the 
main loop. The main loop is entered with current values of A Q> AA Q) 

X is AX i? and the direction of search in the A coordinates (as denoted 
by the variable SGN) and control immediately falls through to the A- 
loop. In the search mode, the A- loop is repeated at most twice for 
each entry from the main loop. The A- loop increments the current A Q 
by AA 0 in the proper direction (A=A 0 +SGN& AA 0 ) , initializes the state 
vector and other data associated with integration of the state equations 
(Section 2.1.3) and enters the integration loop. 

The integration loop uses the Runge-Kutta-Gill numerical integra- 
tion algorithm to integrate the state equations, saving the r,X coor- 
dinates of the free surface trajectory obtained. The integration loop 
exits on any one of the following conditions: 

a) the trajectory passes through the maximum radius of the con- 
tainer, r>R max ; 

b) the trajectory passes through the axis of symmetry of the 
container, r < 0; 
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c) the trajectory passes through itself, dr/ds>0 and dX/ds<0, 

or dr/ds < 0 and dX/ds>0, where s is the arc length coordinate. 


These conditions are monitored by the Boolean function GOBACK and the 
loop is exited with the curfent trajectory in the array SOLN. 


Upon exit from the integration loop, initialization is performed 

for the X-loop and the X-loop is entered. The X-loop is repeated the 

number of times specified by the input variable NX for each entry from 

the A- loop. It*s function is to search each in the range X+^(NX)*AX 

< X-j_< X-i>(NX)*AXi for an approximate value of the error function, ^ , 

where f = [(V u -V Uj . J )/100.]2+ [(0 -9 . )/l80.] 2 , 

r u current u desired J L current desired 

less than the current value. If a smaller error is found, information 

identifying the current solution is saved and the loop is repeated until 

completion. Upon completion of the X-loop, one of two conditions exist I 

either a better solution was found, or a better solution was not found. 

If a better solution was not found, control passes to the top of the A- 

loop and the process is repeated for A=A 0 -SGN*AA 0 , If no better solution 

is found on the second pass, the presumption is that AA Q and AX^ are too 

large; they are halved and control passes to the top of the main loop for 

a new search about the same Aq and Xi as before but with smaller values 

of AA 0 and AX^. 


If a better solution is found upon completion of the X-loop, an 
accurate value of is computed for the current solution, and, if ^ is 
greater than the input tolerance, EPSC, A Q is set equal to A and control 
passes to the main loop for another search with no change in direction. 
This procedure is repeated until either convergence is established 
(tff < EPSC), in which case final output is generated, or until AA 0 be- 
comes less than a prescribed tolerance (currently 10“ 5 times the initial 
input AA 0 ) in which case the message HALVING LOOP, EXECUTION TERMINATED 
is printed and control passes to the input phase for the next case, if 
any. 


Final output consists of tabulated values of the free surface equili- 
brium shape coordinates R and X as well as values of R* , X*, and incre- 
mental ullage volume V* . 


2. 1.6. 3 Container Definition - In order to provide the maxi- 
mum flexibility in specifying a container shape, it was decided to impose 
on the user the burden of writing a FORTRAN subroutine called CAN, with 
entry points RCAN and RCANP. RCAN is passed a value of X as an argument 
and computes the corresponding value of r. RCANP is passed a value of X 
as an argument and computes the corresponding value of dr/dX. The version 
of CAN presented herein (Section 6.1.1) computes these values for the con- 
tainer selected for the vibration analysis demonstration problem and can be 
used as a model for other containers. 
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2.2 Mass and Stiffness Matrices 

Mass and stiffness matrices of the complete structure (only fluid 
is used in this study because the tank wall is assumed rigid) are cal- 
culated using a finite -element approach. In this approach, a continu- 
ous structure is assumed to be composed of simple, small structural 
elements such as tetrahedrons, pentahedrons, and hexahedrons for the 
volumetric fluid elements and triangles and quadrilaterals for the 
surface elements (both gravitational and surface tension). The deri- 
vation to obtain the finite element mass and stiffness matrices is 
based on kinetic energy and strain energy principles, respectively. 

The kinetic energy for a complete structure may be expressed as 


T = h 



P(X, Y , Z) S 


(X, Y, Z, t) dX dY dZ 


(18) 


where T = kinetic energy 


P = mass density 

O 

3 = time rate change of deflection 


t = time 


X, Y, Z = global coordinates 

The difficulty in integrating equation (18) is expressing the deflection 
S (X, Y, Z, t) as a continuous function over the complete structure. In 
the finite- element approach, however, this apparent difficulty is circum- 
vented by idealizing the structure to be comprised of many small structural 
elements for which B (X, Y, Z, t) can be expressed as a continuous function 
within the element boundaries. Thus, the expression (18) is valid for each 
of the finite-elements of the structure. Then the kinetic energy of the 
structure is the summation of the kinetic energies of each of the finite 
elements, that is, 


X = 



T. 

i 


.(19) 


where i refers to one particular finite element n i n . 

The common junction of finite elements is denoted as panel points, nodes 
or joints. Joints will be used here. The deflection S (X, Y, Z, t) is 
easily expressed as a simple function of the joint deflections. These 
element joint deflections are then generalized coordinates or degrees of 
freedom of the complete structure. 
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The approach Is to derive the mass matrix for finite-element, "i", 
in a convenient local coordinate system and then transform it to the 
global coordinate system. The technique is outlined here: 


. (20) 

the mass matrix in the local coordinate system for the 
ith element. This mass matrix is obtained by integra- 
tion using an assumed displacement function. The dis- 
cussion is deferred till later. 

ss the time rate of change of the joint deflections of 
finite-element, "i". This is in the local system. 

The deflections in the local coordinate systems are related to deflections 
in the global coordinate directions by a transformation matrix, [ Y ] of 
direction cosines. Thus, 


where 


[“L ] . 


K (t) } i 


{ h L^oJ’ i bl t i 


( 21 ) 


= the joint deflections of finite element, "i", in the 
global coordinate 

Using equation (21) in equation (20) 

T i = ^ | i M ± 

where [m ] = M? C^] 

\ U 1 i 

is the mass matrix with respect to the global coordinate system for the 
ith finite-element. Further, all the elemental mass matrices are finally 
assembled to give the mass matrix of the total structure, as shown in 
equation (19) . 

The development of the finite-element stiffness matrices is similar 
to that of the mass matrices. The strain energy for the structure may 
be expressed as the summation of the strain energies of each finite 
elements. That is, 


system. 


( 22 ) 


(23) 


where 


{ h G (t) } i 



( 24 ) 
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As was done for the finite- element mass matrix, the stiffness matrix for 
finite-element, M i n , is derived in a convenient local coordinate system. 
Thus , 


U. 

1 


M'Wi M 1 {v^}. 


(25) 


where 



= the stiffness matrix with respect to local coordinate 
directions for finite element, M i". This stiffness 
matrix is obtained by integration, using an assumed 
displacement function. This will be discussed later. 



the joint deflections of finite element, i, measured 
in local coordinate system. 


The same transformation matrix, Mi, which was used in equation (22) 
is used, here to relate the deflections in local coordinates to deflections 
in global coordinates. Substitute then to give 


U i = ^ { h G (t) } i M ± t (26 ^ 

where [kJ ^ = [v]* [\] . [y]_ (27) 

is the stiffness matrix with respect to the global coordinate system for 
the ith element. . 


Euler angle rotations at some joints (where the body coordinate is 
neec j e( j to be different than that of the global coordinates) are input in 
the program to allow the joint degree of freedom at these joints to be 
different than that of global X, Y, Z directions. 


2.2.1 Surface Tension Finite Element - The basic surface tension 
element is a triangle; quadrilateral elements are formed by taking the 
average of the four overlapping triangles created by the diagonals. 

For each triangle element, a local Cartesian coordinate system is 
defined such that vertex 1 is at the origin, vertex 2 is on the positive 
x-axis, vertex 3 is in the positive y direction as shown below. 



y 
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Figure 2-6. Triangle Coordinate System 


The displacement field, in the normal direction (z) , is chosen to 
be. compatible with the triangle gravity element and the tetrahedron fluid 
element and will be linear of the form 


w=a+bx+cy (28) 

The coefficients (a, b, c) are eliminated in terms of the 3 vertex dis- 
placements (w^, W£, W 3 ) . 

The mass matrix for the surface tension element is zero because 
there is no contribution of this element to the kinetic energy of the 
total system. 

The stiffness matrix for the surface tension element is obtained 
as follows. Surface energy is associated with the liquid/vapor inter- 
face and the work done by the external liquid molecules to extend the 
surface. To simplify calculations, a hypothetical tension that acts 
in all directions parallel to the surface is substituted for the sur- 
face energy. This hypothetical tension is generally termed "surface 
tension". Surface tension has the same dimensions as surface energy 
per unit surface area, and it must have the same numerical magnitude# 
The concept of liquid surfaces behaving like a stretched membrane must 
not be misconstrued because surface energy is the fundamental liquid 
property and surface tension is merely a mathematical equivalent. The 
derivation given here is similar to that given in Reference (2) under 
Flexure of Plates with Simultaneous In-Plane Forces. 

The effect of the normal deflection, w, is to introduce additional 
in-plane strains. Considering the figure (2-7) below, if points A and 
B move vertically, then the original length Ax becomes 

|<-> 2 - (s “) 2 | 4 


(29) 
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1 + % 



(30) 



Figure 2-7. In-Plane Strain Due to Normal Deflection 


The extension of the plate is then (neglecting higher order terms) 


w 


u = % AX I ^ 


(31) 


and the in- plane strain (u/ Ax) is 


£ x = * 


■s \ 2 
OV7 

c)x 


(32) 


Similarly, in the y direction. 


_ i i 

e y“Msy 


(33) 


The in- plane strain energy thus becomes 




st 


£) + (t^ 1 y dxdy 


(34) 


Where <j is the coefficient of surface tension. Noting that the slopes 
are related to vertex displacements, we may Vrite 


5w 
dx 
5 w 

Sy 


= [6]{«{ 


(35) 
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where 


M T = [ w l> w 2’ w 3 ] (36) 

and [Gj is a differentiation matrix. 

The potential energy for the element becomes 

VH«r OJld (37) 

where [K st J is the stiffness matrix for the surface tension element and 
is defined as 


[ K st ] - -// [ G 1 T t G ] dxd y 


Using the linear displacement function of Equation (28) 


= [1 x y] b 


To evaluate a, b } c 


0 0 


x 2 0 
X 3 y 3 


From which 
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where 


[A] 


-1 



0 

l/x 2 


^(x 3 -x 2 )/x 2 y 3 -x 3 /x 2 y 3 



w = [X x y] [A] 1 \S\ 


(42) 


(43) 


Now 


r 


c3w 

£x 

Sw 

dy 


0 1 
0 0 


By comparing the above 


0l W l M*l 


equation with Equation (35) , we see that 


(44) 


[<0 


_1_ 

2A 



0 


x 3 -x 2 -x 3 x 2 


(45) 


where A = % x 2 y 3 (the area) . 
Then from Equation (38) 


[ K st ] = "■ A M T M 


4A 

A + (x 3' x 2 )2 

~A ' X 3 (x 3' x 2 ) 

x 2 (x 3 -x 2 ) 


'A " x 3 (x 3' x 2 } 

V 2 + X 2 

y 3 3 

" x 2 x 3 


x 2 (x 3 ”X 2 ) 

" X 2 X 3 

x 2 

2 J 


2.2.2 Gravity Finite Elements - The basic gravity element is a 
triangle; quadrilateral elements are formed by taking the average of 
the four overlapping triangles created by the diagonals. This element 
uses a linear displacement field (’Q) that is boundary conformable. 

The mass matrix for the gravity element is zero because there is no 
contribution of this element to the kinetic energy of the total system. 
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The stiffness matrix for the gravity element is obtained by expres- 
sing the gravitational potential energy in terms of the vertex displace- 
ments as 


U 

g 



area 


n) (w * ■§) ds 


. (47) 


where U = gravitational potential energy 

g 

*h = unit outer normal 

e = a unit vector parallel with the gravity vector ~g, but of 
opposite sense, i.e., ’e = -"g/g 

A noteworthy observation can be made with reference to the gravitational 
potential energy expressed in Equation (47), We note that since we have 
a boundary conformable element, the surface integrals such as Equation 
(47) will all cancel each other throughout the interior of the fluid in 
a container, since n on common element boundaries is equal and opposite. 
Thus, the gravitational potential energy will depend only on displacement 
coordinates at the boundary of the entire volume of fluid (the free sur- 
face and the. wetted container wall). 

Notice also that for a rigid tank, w • n is non-zero only at the 
free surface where e = *n; thus, 

U -\p g / (w • n) 2 ds, (48) 

& ‘■'rree 

Surface 


a much more familiar expression than that of Equation (47). 

2.2.3 Fluid Finite Element - The basic fluid element is a tetra- 
hedron; pentahedron elements and hexahedron elements are synthesized, 
simply by placing six and ten overlapping tetrahedrons together, re- 
spectively and averaging the result. The averaging is carried out to 
eliminate the bias, if any. 

For each tetrahedron element, a local Cartesian coordinate system 
is defined so that vertex 1 is the origin, the x-axis includes vertex 2, 
vertex 3 lies in the x-y plane and vertex 4 always has a positive z- 
coordinate (Figure 2-8). 
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Figure 2-8 . Local Coordinate System for Tetrahedron Element 


This element uses a linear displacement field (constant strain) 
that is boundary conformable. The displacement field throughout the 
element is expressed in terms of coordinate locations and appears as 


w (x, y, z, t) = a Q + a 1 x + a 2 y + z (49) 

The coefficients a k (t) , k ■ 0, 1, 2, 3 are eliminated in terms of the 
12 vertex displacements. 

The mass matrix for the fluid elements is obtained by expressing 
the kinetic energy as 


T = \ 


w • w p dv 


(50) 


Vol 

where T - kinetic energy 

■ • » • • 

w = "a + aL x + 1L y + "a- z 

o 1 2 . J 

P - mass density 

This gives rise to a (12x12) mass matrix. 

The stiffness matrix for the fluid element is obtained by expressing 
the volumetric dilatation strain energy in terms of vertex displacement 
coordinates as 


u D = % 


Vol 


K9 2 dv 


(51) 
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where = volumetric dilatation energy 

K = fluid bulk modulus 

0 = volumetric strain 

2 . 3 Vibration Analysis 

Two methods are included in this study for the solution of the 
eigenvalue/vector (i.e., frequencies/mode shapes) problem. The first 
method uses a Jacobi technique and is documented in Reference 3 in the 
description of Subroutine M0DE1. This subroutine is used to calculate 
frequencies /mode shapes for small and intermediate size problems (ap- 
proximately 120 DOF on a computer with 65000 core). No description of 
the method will be given here because the use of this subroutine is 
straightforward . 

For larger size problems, a second method of calculating the fre- 
quencies/mode shapes is included. This is the iterative Rayleigh- Ritz 
method which is described in References (1) and (4) . In this method, 
a large problem is reduced to a smaller problem in a particular fre- 
quency range. Mode shapes are initially assumed and then by the itera- 
tive technique are improved until they converge to the normal vibration 
modes of the structure. Because the iterative Rayleigh- Ritz method is 
not as straightforward as the first method described above for smaller 
size problems, the iterative Rayleigh- Ritz method is briefly described 
here 


For a discrete coordinate model of a structure having n degrees of 
freedom, the equations of motion can be written as 


M |h} + [K] {h[ = 0 (52) 

where jhj = jh(t)} vector of discrete coordinate displacements, 
M = mass matrix 
M - stiffness matrix 

If a solution of the type |h[ ~ }hj e :LVt , implying a simple harmonic 

motion is assumed, equation (52) can be written as 

([K] - - 2 [M]) \h\ = {0j (53) 

Equation (53) is recognized as a matrix eigenvalue problem of order n, 
whose eigenvectors £ 3> ] are the mode shapes and whose eigenvalues 
are the frequencies. A complete sequence of trial vectors 
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{v} l , M 2 , {v| 3 


W n (54) 


which are linearly independent, is assumed. The displacement jhj is 
then expressed ss a linear sum of the first m trial vector, that is, 


M 

= [V] {q| 

(55) 

(nxl) 

(nxm) (mxl) 


Substi tut ion 

of Equation (55) into (53) 

T 

and multiplying by [Vj gives 

([»] • 

■ “ 2 O*]) = I 0 ! 

(56) 

where 

i — i 

£ 

i i 

= [v] T [k] [v] 

(57) 

and 

[M*] 

= [V] T [M] [V] 

(58) 


Equation (56) is a matrix eigenvalue problem of reduced order "m" whose 
eigenvectors are [**] and eigenvalues are • The solution of 

equation (56) has the form 


M = [**] 


(59) 


where j q*[ is the normalized coordinate vector. The eigenvalues, 
f«5 2 J j approximate the first "m" eigenvalues of the original structure. 
The associated eigenvectors [$] of the original structure are obtained 
by substitution of Equation (59) into (55), yielding 


or 


{ 


M 

= [V] 

[**] h*i 

(nxl) 

(nxm) 

(mxffl) (mxl) 

h} = 

MK1 



( 60 ) 
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where 


M = [v] [y*] 


(61) 


2 

The accuracy of the mode shapes £$1 and frequencies obtained 

depends entirely upon the trial vector . If [V] contains the true 

modal patterns, then the eigensolution for[4>J and[^<i>^] are exact. 
However, in general, that is not the case. Exact results can be obtained 
for the first "m" modes of the structure if the trial vectors [V] do not 
have any contribution from modes higher than n m". Thus, an improved set 
of trial vectors can be calculated by suppressing the contribution of 
higher modes in approximate mode shapes. The procedure for suppressing 
the contribution of the higher modes is well known; in fact, it is the 
basis of the Power or Stodola- Vianello 1 matrix iteration method of modal 
analysis. Here, however, the method is applied to all modes simultane- 
ously and is given as. 


[K] [V] = [M][>] 


(62) 


The solution is carried out for [v] , which is then used to repeat equa- 
tions (56) through (62) . The cycle can be repeated until all the mode 
shapes and frequencies have converged to within a prescribed 

tolerance. Convergence is assured because the technique is equivalent 
to a power iteration applied simultaneously to all modes. Thus, the 
convergence theorems associated with the power method are directly ap- 
plicable. The role of the eigensolution (equation (56)) is to prevent 
all modes from converging on the lowest mode. 

Associated with the iterative Rayleigh-Ritz technique are para- 
meters that affect the convergence and, hence, computer time which will 
be briefly discussed here. They are: 

a. the initial mode shapes assumed to start the iteration process, 

b. the number of modes used, 

c. the repression of higher modes, and 

d. shifting. 

a. Initially Assumed Mode Sha p es - The choice of initial mode 
shapes plays a very important role in the success of the technique. 
Inherent with the initial mode shape selection are two basic problems: 
(1) modes may be missed, and (2) the triple product M - [V] T [h] [V] 
may be ill-conditioned if the columns of [V] are not sufficiently inde- 
pendent. It does not appear that there is a way to guarantee that the 
above two conditions will be met with any selection of [VJ , however, 
the chance of them occurring can be minimized with some judicious se- 
lection of the vectors. If the elements of the vector or of matrix tv] 
are randomly generated, it has been found that the chances of the above 
two conditions being violated is very remote. 
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b. Number of Modes Used - An increase in the number of modes used 
will, in general, decrease the number of iterations required for conver- 
gence. However, if more modes are used, the computer time for each 
iteration will increase because of the increase in sizes of the matrices 
used. Determination of the optimum number of modes to use requires an 
empirical assessment. 

c . Repression of Higher Modes - As pointed out earlier, exact 
results can be obtained for the first "in" modes of the structure if the 
trial vectors in [V] do not contain any contribution from modes higher 
than "m". Generalizing, it can be said that an improved set of trial 
vectors can be calculated by suppressing the contribution from the higher 
modes in the approximate mode shapes at each step. This is achieved as 
follows , 


[v] . = [K] - 1 [m] [V] ,. x 


(63) 


The subscript j denotes the iteration number. If this iteration is 
repeated sufficient number of times, modes corresponding to the lowest 
frequency will be reached. If this iteration is repeated too many times, 
the mode will repeat itself in one or more columns of [V] and will render 
[V] I [H] [V] to be ill-conditioned. The use here is not to converge 

to a mode but just to repress the higher modes and, hence, just a one time 
application is advisable. 

d. Shiftin g - Shifting is an useful technique to speed the conver- 
gence of modes whose eigenvalues are close to the shift value. An addi- 
tional benefit of shifting process is the conversion of the stiffness 
matrix (in case of a free-free structure) from singular to a non-singular 
matrix. The method is as follows. 

To introduce the shift value, A s , the quantity A s [m] is added and 


hf = { 0 } . (64) 


Define 

jt] = [K] - A g [M] (65) 

and 


subtracted in Equation (53) to give 


/ 


[k] - A. [m] - a. 7 [M] + A s [M] j • 


( 66 ) 
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to give 


[£] - ° 2 M ) {hf = {0} 


(67) 


This is 
is 


now the eigen- problem to be solved rather than (53) . 
non-singular even if l>] was not. 


Note that 


The eigenvalues of the original system are easily obtained as 


2 n 2 , . 

G> - i/ + A. 

S 


( 68 ) 


2 

The convergence will be to the lowest absolute value of f2 . Thus, 
shifting by a value, A. s , the eigenvalues, 3 around this shift point 
are converged to first. 

Some general remarks on shifting follows. 

1. Analysis of a Free-Structure - Because a free structure has a 
singular stiffness matrix, the solution of the simultaneous 
equations in the iteration loop is not possible. However, the 
shift technique alleviates the problem. 

2. Specific Frequency Range - When a shift value is used, the 
modes with eigenvalues closest to the shift value will con- 
verge first, which enables one to obtain the modes in the 
desired frequency range only. 

3. Large number of modes - By repeated use of different shift 
values, any number of modes can be obtained. 

4. The following observations are made based on the results of 
Reference (4) . 

7 9 7 

a. If the lowest eigenvalues in the range g>j-, 
are needed, a shift value of zero should 

be used for a restrained structure and one for a free-free 
structure. 

b. If the modes are needed in an intermediate range, a shift 
midway between the lowest and the highest expected eigen- 
values should be used. 

2.4 Finite Element Model 


After the static free surface shape has been established, the total 
fluid is modeled as an assemblage of finite elements for the vibration 
analyses. Instead of analyzing the complete tank and fluid, it is con- 
venient to reduce the number of degrees of freedom (and, thus, reduce 
the computer time) by using a 90° sector. This 90° model requires four 
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sets of boundary conditions to completely represent the total 360° model. 
This technique will be discussed in Section 2.4.1. To reduce the amount 
of input data describing joint coordinate locations, degrees of freedom, 
and Euler angles along with the finite element joint numbers, a subrou- 
tine to generate this data was developed and is described in Section 2.4.2 

The vibration analysis is for an axi- symmetric tank with rigid walls 
as shown in Figure 2-9. Euler angle rotations are used to give body coor- 
dinate systems normal to the tank walls and, thus, allow for fluid slip- 
page tangentially and zero penetration normally. For consistency, the 
V DOF is normal to the tank wall. 


X 



2.4.1 90° Model Boundary Conditions - Using the four sets of 

boundary conditions given in Figure 2-10, the complete 360° model is 
represented. Because this study is only concerned with the lateral 
"y" slosh modes, the symmetric/anti- symmetric boundary condition was 
the only condition used for the data generator. 



Type of 

: Motion 

Boundary 1 
(XY plane) 
Constraints 

Boundary 2 
(XZ plane) 
Constraints 

Center * 

Constraints 

Boundary 1 
(XY plane) 

Boundary 2 
(XZ plane) 

Allow 

Fix 

Allow 

Fix 

Allow 

Fix 

Symmetric 

1 

Symmetric 

8 , 8 

y 

8 z 

V S z 

8 y 

S x 

8 , 8 
y z 

Symmetric 

Anti- Symmetric 

8 x> 8 y 


8 y 

»*■ 8 z 

8 y 

S , S 

X* z 

Anti- Symmetric 

Symmetric 

s z 

V 8 y 

V 8 z 

8 y 

S z 

V S y 

Anti- Symmetric 

An t i - Symme t r ic 

8 

2 

V 8 y 


8 , K 

x z 


V S y’ 8 z 


Figure 2-10. 90° Model Boundary Conditions 
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2.4.2 Data Generation - A computer subroutine (LBDGEN) was developed 
to generate joint coordinate locations, degrees of freedom, and Euler 
angles along with the finite element joint numbers. Using a drawing of 
the container wall and fluid surface in the XY plane, the analyst sketches 
the desired grid and obtains the joint X, Y coordinates and Euler angles 
(at the container wall). With this information, the number of sectors 
desired in 90°, and other information detailed in Section 6.2.2, the data 
is generated to calculate the finite element mass and stiffness matrices 
for the fluid compressibility, gravity, and surface tension. 


To avoid any ambiguity of the mode number of the first slosh mode, 
an algorithm to calculate the first slosh mode number was obtained as 
follows. The total number of degrees of freedom is given as 


NDOF - M + M + M (69) 

o s c 

where M = number of circulation (zero frequency) modes, 
o 

M = number of slosh modes which is the number of surface degrees 
of freedom, 

M = number of crunch (high frequency) modes which is the number 
c of fluid elements. 

Because NDOF, M s and M<, are easily calculated from the finite element 
geometry, the number of circulation modes is obtained as 

M - NDOF - M - M 
o sc 

The mode number of the first slosh mode is then M q + 1. 

In terms of the grid sketched in the XY plane; NDOF, and 
are given as 

NDOF = NGPAX + (2*NGPCW+3*NGPI S ) *N SECT 
M = 1 + (3*NGPFS+2)*NSECT 

S 

M = NGPEL*NSECT 
c 

where NGPAX = number of grid points on X-axis, 

= number of grid points on container wall (except at 
X-axis) 


NGPCW 
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NGPIS - number of grid points on interior and surface (except 
at X-axis and container wall) 

NSECT = number of sectors in 90° model 

NGPFS = number of grid points on fluid surface (except X-axis 
and container wall) 

NGPEL = number of grid point elements, i.e., number of elements 
on XY plane 

Usage of these algorithms is demonstrated in Section 3.2. 
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3 . RESULTS 


3 . 1 Free Surface Static Equilibrium Sha pe 

The digital computer program for generating the free surface static 
equilibrium shape has been verified by running several containers in both 
the survey mode and the search mode. The survey mode results consist of 
tabular printout and several plot frames. The tabular printout, a sample 
of which is shown in Figure 3-1, consists of values of ullage volume and 
contact angle vs tank axis intercept (X) for each value of A specified 
in the input data (Section 6.2.1). The asterisks indicate that no solu- 
tion exists for this value of A and X. The plot output graphically dis- 
plays this data in a manner intended to allow the user to visually bracket 
the desired values of ullage volume and contact angle and thus provide 
good initial values as input to the search mode. Ullage volume percent 
vs tank axis intercept plots are shown in Figure 3“2a; Figure 3- 2b gives 
the correspondence between the numbers associated with each curve and 
the value of A. Contact angle vs tank axis intercept plots are shown 
in Figure 3- 3a; Figure 3-3b gives the correspondence between the numbers 
associated with each curve and the value of A. 

The range on tank axis intercept is XUP- AX to XLCH* AX where AX is 
given by (XUP-XLO)/NX and XUP, XLO and NX are input. by the user. . Figures 
3-4a and 3-4b show the cross plot of ullage volume percentage vs contact 
angle for two. consecutive values of A. The range of A is from ACOFO+DACOF 
to ACOFO+NA* DACOF where ACOFO, DACOF, and NA are input by the user. Figure 
3- 4c gives the correspondence between the numbers associated with each 
curve and the value of tank axis intercept. All of the, above figures 
(3-L through 3-4) apply to the tug like tank shown in Figure 3-8. 

As an example of how to use the survey mode plotted output, consider 
determining initial conditions to the search mode to find the equilibrium 
shape for an ullage volume of 507* and a contact angle of 45°. For this 
particular tank, the cross plots, Figures 3-4a and 3-4b, are the most 
useful. The idea is to bracket the desired point and the solutions for 
A=0. 05 (Figure 3-4a) lie to the left of the desired point and the solu- 
tions for A=0 . 04 (Figure 3-4b) lie to the right. Therefore, reasonable 
input values might be ACOFO-O.04 and DACOF=0.005. To determine the range 
on tank axis intercept note that the desired point lies between points 
numbered 17 and 19 on the plots. Referring to Figure 3-4c, point 17 
corresponds to X=84.2 and point 19 corresponds to X=76.9; therefore, 
reasonable input values might be XUP=84, XLO=77, and NX-20. 

The search mode has been run with a number of containers. Bond 
numbers, contact angles and ullage volumes.. Results are plotted for 
a cylindrical tank (Figure 3-5), a spherical tank (Figure 3-6), a cosine 
tank (Fi-gplfi? 3- 7 ) and the Tug- like tank (Figure 3-8) selected for vibra- 
tion- analysis . These results are presented to demonstrate the wide range 
of axis ymme trie- tanks to which the program may be applied. 
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Figure 3-1. SAMPLE TABULAR OUTPUT FROM FREE SURFACE 
EQUILIBRIUM SHAPE PROGRAM, SURVEY MODE 


.AGO 
PHI TO 


1 4 4 .U 
128. 2 
I1S.5 
10 4. 3 
94 .0 

8 M. 0 
7 * .5 

65. 0 
55 *E 
4 6. 4 
36 .9 

2 7.2 
16 .7 
1 3e 2 

3 3 .2 
I 3. 2 
I 3 .2 
1 3. 2 
13 .2 
I 3. 2 
13 «2 
I 3. 2 
13 .2 
1 3. 2 
13 .2 
I 3. 2 
13 .2 
1 3. 2 

7 .7 

12. n 

5 .0 
19. 7 
27 .6 

3 3. 6 
39 .0 

4 3. 8 
48 .3 

5 2. I 

S6.2 



35 


ikHiiliiia 
iibiiiibbiH 
iiibiiibibiibim 
iiiiiimimiiia 




■iliiiiiiiiiiiiiiiiiiiiiiiiiiaiiBillllllll!iii!!8!!!!!!8 


'iiBBiiiiiiiiBiBmi^^lMBflflBBBBBiBmmmBimflflBiflBBBiiBm 

maBBasB^^isaaaaaiiiaBiM 


bbbbbbbibibbib! 


IIBBBBBBB 


BBBIIIBflll 




■BIB! 

BBflfli 


BBBBBBBlBBBBBBBBlBBlBBHiiilii*Bl««i**>**555* J2*55SiJjiS3S*S25S55a555i5SS5S 


flBBBBBBBBB.BBBBBfll 


MMamaa^MBMaaaBaiiaBaBBBWBBMBM 


130 140 


Aais INTERCEPT 


Figure 3- 2a. ULLAGE VOLUME VS TANK 
AXIS INTERCEPT 


p LfH SYH80LS FOR 

ftGE VOtiWff PERCENT ^5 TANK SKIS INTERCEPT 

SYHBOt ' ’ ' S' COP” 


6 

t 

3 

n 

5 

' s 
? 

9 

3 

§0 
& £ 

8 2' 
B 3 
fi *a 


o B QOQOOQP < 00 
". 90000000- Oil 
0 39993999 -01 
«6 9 9 99 9 9 9- 01 

0 soooaooo -in 
0 ^nnuodaQ~os 
O «jononooa -in 

o' 30000000- Of 

o 200a00DU“01 

0 93399999- 02 
e QOOOOaDQ 

- s 99 3999'99-02 
-o 20GQQOOQ -Q l 

- osoaoaooo-'ni 


Figure 3-2b^ ULLAGE VOLUME VS TANK AXIS 

lteRCEPT'7 plot‘"symbols 


3 


liiliiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiigigggggigggiiliggggggggggggggggggigg 

IIIIIIIIIIIIIIIIIIIIIBIBIIIIIIIIIIIlllgllllllllllllllggllllllllllllllglgl 

Illllllllllllllllllllllllllllllggllggllggglglggglglgglllggggglgggllggl^l; 

iiiiiiiiiiiiiiiiiiiimiiiiiiiiiiiiiiiiiigiiggggiigggigigggggggiiigggyii 

lllllllllllllllllllllllllllllllllllllllllllllllllllllllll!il!l!!!!!!!!^ll 


llllllllllllllllllllllllllllllllllllllllllllllllllllllggiiggigigggigl 
iiiiiiiiiiiiiiiiiiiiiiiiiiiiiaigggiiggggggggggggggggggggiggggggggggil 
liiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiiimmiRifl 

!!!!i!!!!!!i!i!!!!!!!!!i i!i!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!linjjl 

■!!!!!!!llil!!i!iii!!!!lii !! !!mj i!i iiii ii i!!!i!!lliiiiliiiii»M 


I 

pi 

II 
fill 
■in 
«ni 
inn 
mil 


llllllllllllllllllllllllllllliri f m2lllllilllllllllll!!!!!l!i!!&M/iil!!l 


III 

mi 


iiiiiiiiiiiiiiiiiiiiiiniiiiiw niiiiiii 
^■HfinannHiHiiiiiiii 


nini!s;n:i 


iniimiii 
punimii 

isiiaaiRaBiiiiaiiiaiiiiiigggggggggg^Qaagggggggggg 
■■■■■■■iiiiiniaiiiiiiiiiiEiiiiiiiiiiiiiiipnniiiiuiiii 

iiiiivaiiriRaBiiiiiBiiBiiiBigiiggiiiggiir^^igiiiiiggigg 

itaniis»iiiiiiisin![/iniiiiiiniiiiiiiiiiiiiiiiiii|nngiiiHiii|ii 

iii^sniiisiiiiininiiRnaiiiiiiiiiiiiiiiiiiiiiiiiiiKinniiiiiiiiiiiii 

llll^s$!i:iilli:!lllll^iilllllllllllllllllllllll||||l^l!ll!l!l!!llll!l 


■52581 


lll|IIIIRai!lilKaiHIKSEIIIIIIIIgllllllllllggggggggl9Sgg|gggggg!ggggg;gggggi 

laiiiiiiiiiiKOKiiiiiiiiiiiiaiBiiiiiaiiiiiiaiiiiafaQBBggggaaiiBgBiggiiggag 


1 1 1 1 1 1 1 III fl 1 1 1 IBI lillBi 


ifilEIIIIIIIIIIIIIIIIIIIIIII 


imiiminiimiiiiiiimmiiiiiiiiniimiBiiiiiiiimiiiiiiiiiiiiiiil 


M 100 >10 120 130 HO 


TANK AXIS INTERCEPT 


Figure 3- 3a. CONTACT ANGLE VS TANK AXIS INTERCEPT 


Q ^ yg ^ gv W* ^ 


38 - 


PUn SYMBOLS FOR 

"COWI ACT AW61 E Of GR IF ES VS 7 $NK AS IS INTERCEPT 

S VWS GL 6 C 0 F 


8 o 60000000*00 

2 .901X00000-01 

0 ?9933399 -01 
. 69999999-01 
.60000000 -06 
. 50000000-01 
. 4 OODOQQO -0 6 
. 30000000-01 
. 20000000 -0.1. 
.99999999-0 2 


M .. . .00000000 

Ti -. 99939999-02 

S3- 2DQOOOOO -0 I 

m -.30000000-01 


Figure 3- 3b. CONTACT ANGLE VS TANK AXIS INTERCEPT, 

PLOT SYMBOLS 



39 



Figure 3-4a. ULLAGE VOLUME VS CONTACT ANGLE , 

A = 0.05 
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Figure 3-6. FREE SURFACE STATIC EQUILIBRIUM SHAPE, 
SPHERICAL CONTAINER 
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3 . 2 Vibration Analysis 

The Tug- like container of Figure 3-8, with a fluid of nitrogen 
tetroxide (N 2 O 4 ) , was used for the vibration analysis. 

From Table A-l of Reference (5), the density (p) and surface tension 
(a) of this fluid are obtained as p = 1.454 dyne sec 2 /cnA ( 1 .. 36x10-4 lb- 
sec2/in4) and <r = 27.4 dyne/cm (1.56xl0" 4 lb/in). From this same refer- 
ence, the contact angle is 0 °- 2 ° for the liquid, its vapor and titanium. 
From Martin test work, the bulk modulus is estimated at 90,300 N/cm 
(131,000 psi) . 

The acceleration (g) is calculated from the equation for Bond num- 
ber which is 


B 

o 


P § r max 


Using the tank radius (r max ) as 81.28 cm (32 in) and a Bond number of 1 
gives g * .00285 cm/sec 2 (,00112 in/sec 2 ) which is 2.9x10-6 Q f the gra- 
vitational acceleration at the Earth s surface. 


Note that there is an inconsistency in the characteristic length 
used here (r max ) and that used in the static equilibrium section 2.1.4 
0W>* Because = 372.11 cm (146.5 in), the acceleration used^in 

the vibration analysis should have been reduced by (372.11/81.28) - 

20.96. Because the terms of the gravitational stiffness matrix are al- 
ready a factor of 10"* 2 less than the terms of the surface tension stiff 
ness matrix, further reduction in the gravitational stiffness matrix 
terms should not effect the vibration analysis results that were ob- 
tained . 


Four different size models were used in the vibration analysis ^ 

A course grid with one sector size and a fine grid with three sector 
sizes. The grids are shown in Figures 3-9 and 3-10. The algorithms 
developed infection 2 . 4.2 are used here to calculate the total number 
of degrees of freedom (ND0F) , the number of slosh modes (Mg) , the num- 
ber of "crunch" modes (M^,) , and the number of circulation modes ‘(Mq) 
in terms of the grid point geometry. These results are given in Table 
3-1. The mode number of the first slosh mode is simply M Q + 1. 
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Table 3-1 


NUMBER OF DOF, SLOSH MODES, CRUNCH MODES, AND 
CIRCULATION MODES 



Grid 1 
2 Sectors 

Grid 2 

2 Sectors 

3 Sectors 

4 Sectors 

NGPAX^ 

4 

5 

5 

5 

NGPCfrK 1 ) 

6 

8 

8 

8 

NGPIS^ 12 

4 

11 

11 

11 

NGPFS^ 1 ) 

4 

5 

5 

5 

ngpel(D 

7 

15 

15 

15 

nsect(D 

2 

2 

3 

4 

— — — — 

— --- — — 

• - - 

- ; — — — 

’ ' • 

ND0F 

52 

103 

152 

201 

M s 

29 

35 

52 

69 

Me 

14 

30 

45 

60 

Mo 

9 

38 

55 

72 


(1) Symbols are defined in Section 2.4.2 


Using nominal values of bulk modulus, acceleration, fluid density 
and surface tension, a stiffness matrix composed of 

[ K fluid (BKM) ] + [ K gravity (s ’ P } ] + [ K SURFTN (tr) ] 

with terms on the order of 10^, 10“^, 10“ 5 respectively with Grid 1 and 
2 sectors in the 90° model. Based on frequency, there was no problem 
separating the high frequency "crunch" modes. However, it was difficult 
(based only on frequency) to identify which frequencies were slosh modes 
and which frequencies were circulation modes (should be zero) . Factors 
of 10 2 , 10^, 10° were applied simultaneously to acceleration and surface 
tension. A summary of results (Table 3-2) shows the to 2 of the slosh 
modes to vary directly with the factor used and the o 2 of the circula- 
tion and crunch modes to be uneffected. Scale factors were also applied 
to bulk modulus. A summary of results (Table 3-3) shows the tu 2 of the 
slosh modes to be uneffected but the w 2 of the circulation modes de- 
creased almost directly with the factor used and the w 2 of the crunch 
modes decreased directly with the factor used. 
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Table 3-2 


EFFECT OF SIMULTANEOUS ACCELERATION AND SURFACE TENSION 
VARIATION ON FREQUENCY {<£) 

(E-7 = 10“ 7 ) GRID 1, 2 SECTORS 


Mode Type 

Mode No. 

Nom g,<r 

NomxlO^ 

NomxlO 4 

NomxlO^ 

Circulation 

9 

2.578E-7 

1. 93 IE- 7 

6.335E-7 

3.638E-7 

Slosh 

10 

1.403E-6 

1.289E-4 

1.289E-2 

1.289E-0 

11 

2.499E-6 

2.389E-4 

2.390E-2 

2.390E-0 


37 

9.765E-2 

9.765E-0 

9.765E+2 

9.763E+4 


38 

1.194E- 1 

1.194E+1 

1.194E+3 

1.193E+5 

Crunch 

39 

1.453E+7 

1.453E+7 

1.453E+7 

1.454E+7 


Table 3-3 


EFFECT OF BULK MODULUS VARIATION ON FREQUENCY O 2 ) 
(E-7 = 10“ 7 ) GRID 1, 2 SECTORS 


Mode Type 

Mode No . 

Nom BKM 

NomxlO ^ 

NomxlO" 4 

NomxlO -8 

NomxlO" 8 


Circulation 

9 

2.578E-7 

3.299E-9 

1 . 118E-10 

2.633E-13 

2.590E-15 

5.307E-17 

Slosh 

■ 

1.403E-6 

1.293E-6 

1.289E-6 

1.289E-6 

1.289E-6 

1.285E-6 


mm 

2.499E-6 

2.392E-6 

2.390E-6 

2.390E-6 

2.390E-6 

2.370E-6 


■ I 

9.765E-2 

9.765E-2 

9.765E-2 

9.763E-2 

9.565E-2 

1.023E-2 



1 „ 1.94E- 1 

1.194E-1 

1.194E-1 

1.193E-1 

1.148E-1 

1.297E-2 

Crunch 

39 

1.453E+7 

_i 

1.453E+5 

1.453E+3 

1.454E+1 
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Scale factors were also applied to bulk modulus for Grid 2. A summary 
of results (Table 3-4) shows the same variation of a?- with scale factor 
as was obtained with Grid 1. 


Table 3-4 


EFFECT OF BULK MODULUS VARIATION ON FREQUENCY (o>^) 


(E-7 = 10" 7 ) GRID 2, 2 SECTORS 


Mode Type 

Mode No. 

Norn BKM 

NomxlO ^ 

NomxlO ^ 


Circulation 

38 

2.980E-6 

3.603E- 10 

3.403E- 14 


Slosh 

39 

5. 545E-6 

2 . 072B-6 

2.072E-6 

2.060E-6 


40 

6.055E-6 

4.829E-6 

4.829E-6 

4.816E-6 


72 

2.840E-1 

2.840E-1 

2.185E-1 

1.329E-1 . 


73 

4.309E-1 

4.309E-1 

2.324E-1 

1.356E-1 

Crunch 

74 

1.082E+7 

1.082E+3 

2.794E-1 

1.545E-2 


In addition to frequency, the effect on modal displacement with variation 
of bulk modulus was also noted for Grid 2 with the results summarized in 
Table 3-5. 


Table 3-5 


EFFECT OF BULK MODULUS VARIATION ON MODAL 
DISPLACEMENT. GRID 2, 2 SECTORS 
MODE 39 (FIRST SLOSH MODE) 


Grid Point 
Number 

r~ 

Nom BKM 

— 4 

NomxlO 

NomxlO ^ 

NomxlO*" 10 

20 (8X) 
24 (8X) 

-1.193 

5.803 

-.2069 

-2.814 

- .2069 
-2.814 

-.2075 

-2.788 
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Based on the results of the above four tables, it was decided to use a 
scale factor of 10“ ® on bulk modulus to assure valid frequencies and modal 
displacements for the slosh modes. With this scale factor. Table 3-6 shows 
the variation in the first three slosh modes with the various grids used. 
The plotted mode shapes for the first 3 slosh modes are given in Figures 
3-11 through 3- 13 for Grid 2, 3 sectors. The undeformed and deformed _ 
surface joints are shown in a perspective view. Inclusion of the internal 
and wall joints was tried but the large amount of plotted data made the 
viewing too difficult and was thus abandoned in favor of just showing the 
surface. Plots of only the joints in the XY plane was also tried but did 
not give as satisfactory a plot as the surface perspective plots. Capa- 
bility exists in the computer program for any of these plots at the user 
option, however. 


Table 3-6 


SLOSH FREQUENCY (HZ) 
NOM BKM*10~ 8 


Slosh Mode 

Grid 1 
2 Sectors 

Grid 2 

2 Sectors 

3 Sectors 

1 

4 Sectors 

1 

2 

3 

.0001807 

.0002460 

.0004080 

.0002291 

.0003497 

.0006205 

.0001841 

.0002404 

.0003224 

.0001600 

.0001605 

.0002299 

L. — — — 




LEFT EVE VIEW 
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MODE 4, F = .000194 HZ. 

CENTER OF EYES LOCATION 
X = 7. OOOOOOOOE+Ol 
Y = 0 . 

Z =-2. OOOOOOOOE^Ot? 

RUN NO. = G2/3-S 


I US OX TANK. 9N=1. ULL VOL = 00. GRID 2, 3 SECTORS 
VIEW POINT LOCATION ROLL ANGLE = 160.0DEG 

X = 5.00000000E+01 CONE ANGLE = 30.0DEG 

* * 0- EYE TO EYE = 5.0 IN 

Z - 0. 

DATE = E4MR75 


Figure 3- IX. FIRST SLOSH MODE, SURFACE PERSPECTIVE VIEW 


LEFT EYE VIEW 
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MODE 5, F - .0002HO 

CENTER OF EVES LOCATION 
X = 7.00000000E+0! 

Y = 0. 

Z =-2.00000000E+02 
PUN NO. = G2/3-S 


. I US OX TANK. BN=1. ULL 
VIEW POINT LOCATION 
X = 5.00000000E+01 
Y = 0. 

2 = 0 . 

DATE = 2HMR75 


VOL = 90. GRID a, 3 SECTORS 

ROLL ANGLE = 180.0DEG 
CONE ANGLE = 30.0DEG 
EYE TO EYE a 5.0 IN 


Figure 3-12. SECOND SLOSH MODE, SURFACE PERSPECTIVE VIEW 



LEFT EYE VIEW 
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MODE G, F = .000322 HZ. 

CENTER OF EYES LOCATION 
X * 7. QQOGOOOQE+O 1 
Y * 0. 

Z =-2. OQOOOOOOE+02 
RUN NO. = GE/3-S 


IUS OX TANK. BN=I. ULL VOL = 90. GRID 2, 3 SECTORS 
VIEW POINT LOCATION ROLL ANGLE = IQO.ODEG 

X = 5.00000000E+0 1 CONE ANGLE = 30.0DEG 

Y * 0. EYE TO EYE = 5.0 IN 

Z = 0. 

DATE = 24MP75 


Figure 3-13. THIRD SLOSH MODE, SURFACE PERSPECTIVE VIEW 
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4. CONCLUSIONS 

It is felt that a significant contribution was made by this study 
to the state of the art in finite element fluid analysis at low Bond 
number. In this study methods and computer programs for definition of 
the free surface static equilibrium shape at low Bond number, calcula- 
tion of the stiffness matrix due to surface tension, generation of 
joint coordinate locations, degree of freedom values, Euler angles, 
and element joint numbers, and calculation of the vibration mode shapes/ 
frequencies and plotting of these mode shapes were derived and coded. - 
As with probably all new investigative analytical studies, review of 
the work performed reveals that a "blind-alley" was investigated and 
that there are several items that should be studied further. 

To determine the free surface static equilibrium shape, an energy 
minimization technique was originally attempted. In this method, the 
displacement state was sought for static equilibrium corresponding to 
the minimum potential caused by gravitational potential energy, surface 
tension potential energy and the virtual work done by ullage pressure 
acting through virtual displacements, all subject to the constraints of 
constant volume and contact angle at the container boundaries. This is 
the method that was outlined in the proposal for this study, Reference 
(7). Unfortunately, no results were obtained by this original method 
because the attempts at solution continually diverged. Thus, this ap- 
proach had to be abandoned in favor of the force balance method des- 
cribed in Section 2.1. 

One item that should be investigated further is the separation of 
the slosh modes from the circulation modes and "crunch" modes. One 
possible approach, as described in Reference (8), is to describe all 
the joint coordinates in terms of the surface and ignorable coordinates 
by means of constraint equations. With this relationship, the original 
mass matrix is reduced to only the surface and ignorable coordinates. 
This reduced system is further reduced to only the surface coordinates 
by expressing the surface and ignorable coordinates in terms of only 
the surface coordinates. By this technique, only the modal properties 
of the surface slosh coordinates are calculated. 

The data generator subroutine, used to calculate joint X, Y, Z 
locations, degree of freedom values , Euler angles, and finite element 
joint numbers was coded for the lateral slosh boundary conditions, 
that is, symmetric/anti- symmetric boundaries as defined in Figure 2-10. 
Expansion of the data generator to include the symme tric/symme trie and 
anti- symmetric /anti- symmetric boundary conditions should be done. 

Investigation into the non axi- symmetric acceleration field should 
be performed in the definition of the free surface static equilibrium 
shape. This will allow complete generality for the low Bond number 


PRECEDING PAGE BLANK NOT ETT.imq 
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problem. A non. a xl- symmetric acceleration field will require a new data 
generator because a 360^ model definition would be required. The gravity 
stiffness matrix currently allows a non a xi- symmetric acceleration. 
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6. COMPUTER PROGRAMS 


The finite element solution for liquid sloshing at low Bond number 
is accomplished in two main steps. The first step is the free surface 
static equilibrium shape definition and the second step is the vibration 
analysis. Computer programs have been coded for these steps and are listed 
in Section 6.1. Input data to the programs are explained in Section 6.2 
using a sample problem listing. 

A schematic flow chart of the analysis steps is given in Figure 6-1, 
and a brief summary of the important subroutine functions are presented 
in the following pages. 


gRTymni NG PAGE BLANK NOT flLMES 
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Figure 6-1 


FLOW CHART, LOW BOND SLOSH PROGRAMS 
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Summary of Programs : 

STATIC FREE SURFACE obtains the free surface static equilibrium 

shape 

VIBRATION ANALYSIS obtains the vibration characteristics of 

the system (frequencies and mode shapes) 


Summary of Subroutines (used in the VIBRATION ANALYSIS Program) : 

LBDGEN automatic generation of joint X, Y, Z 

values, DOF numbers, Euler angles, and 
element joint numbers for FINELB 

FINELB generates mass and stiffness matrices 

FLUID generates mass and stiffness for fluid 

only 

GRAVTY generates gravity contribution to stiff- 

ness matrix 

SURFTN generates surface tension contribution to 

stiffness matrix 

M.ODED obtains frequencies and mode shapes (small 

size) 

MODES obtains frequencies and mode shapes (large 

size) 


PLOT 


plots the mode shapes 
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6 . 1 Computer Program Listin gs 

Listings of the free surface static equilibrium shape computer 
program and associated computer subroutines is given in Section 6.1.1. 
Listing of the Vibration Analysis Computer program and associated com- 
puter subroutines is given in Section 6.1.2. 
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6.1.1 Listing - Static Free Surface Computer 

Program and Associated Subroutines 


F M3UGGRXN2Q7*TPF* „ LO&ONQ 


IMPLICIT 0 0U3 L £ PRECISION 


PROGRAM 10 DETERMINE HE STATIC E QUTLI 9 R1 UM SHA PE 
FREE SURFACE OF 6 FLJ I 0 WITH A LOW MONO NUM PEP IN 
GIVEN AxlSfMMETRIC CONTAINER 


OF 

A 


HE 


ccccctcccccccccccccccccccccccccccccccccccccccccccccrcccccccccceccccrcecc 

c 
c 
c 
c 
c 


in 

C 



1 1 


LOSTCM. 

GO 4 ACK *DU WHY *A TCE PT «PRI N T *5i A RC 9 

1 2 


RF At 

SURVPL *P HIT 0*P h I CO* PH IF 0* VUP CT S * 0 AS H 

IT 


EQUIVALENCE ISOLNI 500) > »5URVFL 1 1 > ) 

i 4 


01 MENS ION 

SURVPL (51 *5 1 * 2) 

l c 


COMMON 

SOLNt IUU0.5 .2 I 

1 6 


COMMON 

/ CONSTS / P l» ANUM 

IT 


COMMON 

/PARAMS/ A C 0 : *9C0F * 2MA X * RMA X 

1 4 


COMMON 

/ GPRKT •/ QR K (5 I * PR K 1 4 ) 

19 


COMMON 

/TIMES S / DELTA T*T 

20 


COMMON 

/VECTOR/ Y (5 1 *Y OT (61 

2 I 


NAMELIST 

/I NO A 7 6 / AC OF D ** DND NO *D A f. OF *D£ l TA 1 .£ FSC »NA • 

2 2 


i 

IPRNT* NX * PHI 0* PR I NT* RM HU *S E A3 CH» TV 01 • UL P CT* 

23 


2 

XLO *XMA X*XliP 

2 4 

C 



25 


DATA 

NI T/5 /*NOT /5 / *Nt Q/6 / *1 PRNT / HJ / 

25 


OAT A 

tjfl t at /n.iocon/* epsc /i .0 n-trn *ps t si .oo*ns/ 

27 


DATA 

QASH/1H-/ 

29 


OAT A 

PRINT/ .TRUE./ 

?3 


DATA 

A COFO /Q .ODUD/ *D AC OF /O.IDOO/.3CN3NC/1 .HB-OD / *NX / 10/ 

in 


OAT 6 

PHI 0/ 90 ,0 010/ #UL»CT / 50 *0 000/ *X t 0/- 1 .0 010/ 

5 t 


DATA 

RMA X/ -1.0300 / 

3? 


OAT A 

X M «X /- I .0 010/* XUP/- 1 .0000/ 

33 


DATA 

SEARCH/. TRUE./*NA /ID/ 

34 

C 



35 

C 



?fi 

c 

INIT T AL 12 AT TON AN 0 INPUT 

37 

c 



39 

c 



33 

c 

INI T TAL 1 ZE 

CONSTANTS 

40 


ANUM " OHLf 1 AT AN? U .9 • 1 .0) / 45*0) 

4 1 


PI r 1 9D * ODD D * ANUM 

4 ? 


PRKC i 1 “ 3 

. 6 UDO 

4 3 
4 4 

45 

46 

c 

PRK ( ? } z 1 
PRKMI = i 
PR<(4) : 0 

READ INPUT 

.0300 - 3SQRTID.E300J OllTfTXTAT, 

.0 010 * OSQRT (0.5 0301 JjJwiNAX* PAGE IQ 

4 7 

moo CONTINUE 


4 9 


CALL ST ART 


4 3 


RE A D'f NI T * 1 NO A TAJ 

5(1 


N2 = NX 


51 


2 UP - XUP 


5 7 


ZLO = XLO 


53 


2 M A X - XMA 

X 

54 

c 

I NIT I ALU 7 t 

TT ER AT IO~N P AR AMFT £RS 

55 


IFI I N A . S T . 

50 .OR. NZ.GT.50) .AND. . NQT . SEARCH) GO TO 9020 
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56 

5? 

53 

59 

60 
St 
6 ? 
SI 
6 «i 
S 5 
66 
S? 
69 
59 

?n 

71 
7? 
77 
7* 
7 5 

76 

77 
79 
79 
30 
9 t 
9 ? 
9? 
99 

95 

96 
9? 
33 
9? 
90 
9 t 
9? 
33 

9«j 

35 

96 

37 

99 

93 

eon 
1 01 
so? 

103 
1 09 
) 05 
$06 
107 
10 9 
I 09 
I to 
v i t 


NUM6 r ? 

IFt.NOT* SF.4RCH) NUMA = N A 

iF4PM«.LT»n 0 nooo .or . zhax ,lt .o.n nooi go to 

IF 4 ?IJP »L !« 0.0300) ZUP = ZM4 X 

1FI7L0 .LT, 0 .0 0)11) ZLO =0.0000 

04 C OFO = oa C OF 

«COF = MONONO / <7MAX*ZM«I 

7 ^ OL = VUll 13.0000*0.0000 *ZMA X ) 

IflHOWNZ*?) .Nf . 0 1 N 2 = NZ » t 

OZI = I 2UP-ZL0 ) OHt IF IQ4 KNZM 
Z! = (ZUP + ZL 0) f ? .0 DUO 
PSI f r PSI 
SCCFPT = . F SL S t . 

55 M " 1.0000 

C WRITE nftTA 

CALL °A GE HD 
WRIT FIN07. I NO AT At 


C 

c 

c 

c 

c 

Z 003 


c 

c 

c 

c 

c 


31 00 


c 

c 

c 

c ' 

3200 


MAIN LOOP - ITERATE WITH NFW «COFO AND OACOF 


CONT TNUf 
CALL PAGE HQ 
NL OOP = N7 * I 


6 tOOP - GET NEW TRAJECTORY 


no 3000 K A = I • NJ M A 
H r. K A 

T Ft. NOT . 5 f ARCH) i * = \ 

IF j J4 .F 0. ?) SEN = -SEN 
A CCF = ACQFO * S GN ® tl A CO F 

I F | . NO 1 « SEARCH) SUR VPL t) *K4 « ) * ) > = A £ OFT 

IFTPRINU WRITE (NOT * ?0Q 9 ) 4 G AS H* 1=1 e 60S e ACO F* I HAS H* I =t e 50 S 

tin : a.cnou 

V t ?) =0.0 DUO 

7431 = 1.0300 

v <«) = d .u ono 

?!5I : 0.03UD 
nt = n 

7 = O.OQ DO 
00 *100 II - t * N f Q 
ORK I II I = 0*0)00 
CONT TNU f 

'•■ CALL SA VE I N T* I *N£ Q * I A ) 

C ALL TOOT 

IFTPRINT) CALL PR IN TSI NT* ) *NE 0*1 A ) 


i5®®M8 a 


INTEGRATION LOOP* INTEGRATE CI.E.Q. f RON 0 
UNT TL GO HACK = ,F fttSF . 


CALL RUNKT AIN£9*NI» 

CALL SA VEINY* )*NEO*I» ) 
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I 4 2 
I 1 3 

I 1 <i 
I 15 
4 16 
4 17 
I i 3 
4 19 
1 2D 
>2) 

1 22 
123 
1 29 
125 
I 26 
127 
J 2 9 
I ? 9 
t 30 
13 1 
137 
1 3 3 
139 
135 
4 36 
137 
439 
1 39 
MCI 
4 9 4 
IH 2 
1 9 3 
1 M 
M5 
1 9 5 
1 9 7 
t 9 8 

a n 9 
15D 
1 5 < 
15? 
453 
159 
155 
1 5 5 
457 
153 
45 9 
I 5 D 

4 61 
46? 
163 
169 

5 65 
I Bo 
167 


IMPRINT »AN0, MO DINT* IPRNT! • EQ.Ol 
4 CALL PRI NISI NT« 1 *NF. Q ■!* I 

IFC60 8ACK<DUHMM) 60 TO 3200 


C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


92UD 


c 

c 


e 


1 


I 


9 25D 



FINO OF INTEGRATION LOOP 
IFIPRINM CALL PR INI S <NT «4 ® NE3 • Ml 


LOOP ON 21 AT FIXED AC OF TO FINO SMALLER ERROR 
IF ONE FX 1ST 5 

LOOP ON ZI A I FIXED A C OF TO FIND SMALIER ERROR 
TF ONE FX 15 IS 


ZI : 7 ♦ 13LE MLOA 14 NZ/Z ) MOZI 

IFIPRJNT) WRIT F4N0T •2111121 
DC 9H0D I? : l.NLCOP 

Z T = Z I - GZ I 

IF ( e NO T m SEARCH .ANO, KA 0 EQ.») SUR V P L 4 1*2* 1*4*1 I “ ZI 
TFfZI.Gl.ZUP .OR . ZULT.ZLOI 60 TO 9300 
J z a 
CONT TNU £ 

J r J * 1 

I F ( J .G1 • NT 4l ) GO TO 9 300 

IF f RC AN f SOLN'l J* 2 vf A )*ZI I .ST« SOL N 4 J* I *1 A > I SO TO 9200 
HAVE HRACKETEO CAN® DETERMINE PR EL IM 1»N AR¥ VALUES OF 
UL L AG F VOLUME* CONTACT ANGLE* AN 0 ERROR 
VUPC I- VULL tSOLNIJ*5*|A)*SCLN4J*2*IA >4 ZI * ZMA X >• lOU. / TvOL 
P HIF=OHL El AT AN 2 ( S N Gl <$01N< J* 3 . IAl 1*S NGL 4S0LW<J*9* I Al 11 * 
PHIC - 3A TANIRCA NP4 SOLNt J*2 *IA MZI l> 

P HIT - PHJF - PHIC 

IF(FHI= ,L1. 0.03QU • A NO • PHIC .6T « U.OOUUl 
P H II = P HIT 4 2 «0*P I 
PHITD : PHI T /AN JM 

P HI CO = PHI C/ ANU M 
PHIFJ : P H I £ /ANUM 

PSIN - < (VUPCI-ULPCT l/tUO.UOOO) **2 

I (PHI TO -PH 1 3 W 19 0 a 00 00 1 *»*? 

VUPCT S - SNGL 1VUP CT1 

IF | PRINT » UR I TE f NO T.20U7) ZI *P SI N# VU PC TS *PHf T3 * 

PHIfCoPHICC 


IF I SEARCH) 30 TO 9260 
5 UR V PL 1 TZ *1 *K A4* * 11 : 

SURVPL 1 1 2* t *KA« l o2 } z 
GO TO 9000 
CONTI NUE 

T FIPS IN » GT * PS IT l GO 
SAVE VALUES INDICATING 


Z NEW 
PSI T r 
JS AVf 
ISA VE 
AS AVF 
SS3N 
GO TO 


- 2 I 
PSI N 
= J 
r I A 
= ACOF 
: S3 N 
90Q0 


VUP CT 
PHI TD 


TO 9000 
C U R RF. N 7 


MINI WJM E RR OR 
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&SS 9300 CONTINUE 

9 G ^ C NO S GLUT ?ON FOR THIS lACOF.ZI) 

HO IF (PRINT? WRI IE 1 WO T *2001 ) 21 

171 T F 1 S f AR CHI GO TO 90U0 

17? SURVPL II ZH*Ut » .*> = -10110,0 

ui S URVPl < TZ *1 cK A.I f 2) - - J 000 .0 

17 9 A ODD CONTINUE 

8 75 C 

ns c 

177 C FNO OF 7 LOOP 

8 7 9 C 

a 7 9 c 

J9EJ I F 8 S£ 8 RC H 1 SO TO 33110 

ISS ACOFD - flCOF 

192 IMPRINT) CALL PA SEND 

88? go to man 

IBA 3300 CONTINUE 

185 1 F< PS IT *L T ,PS I .AND* ACCEPT) GO 10 2203 

) S S I F 8 I & ,F 0 ■ ?> WR .T TE t NO T *201 I ) A C OF 0 .3 4 C 0F 

187 I F t I 4.EQ ,1 .AND. ACCEPT) GO TO 2)03 

t 88 ACCEPT = p TRUE « 

189 3QO0 CONTINUE 

190 IF1.N0T. SEARCH) 30 TO »50U 

3 98 C 

I 3? C 

193 C F N n OF A LOOP 

1 39 C 

i 95 C 

)3S ? 10D CONTINUE 

197 IF 1 nAHSt 0 ACOF ) .IT. 0 ACO F0.4 .0 0 05) GO T 0 9000 

B93 3 ft C 0 F r 09C0P/2.03U0 

8 99 DZ T = 02 1 / 2.0 non 

200 2 = ZNE M 

201 SGN r -SGN 

?Q2 ACCEPT r . F ft L SE * 

203 GO TO 2000 

209 C GET MORE ACCURATE* VALUES OF ERROR. ULLAGE 

205 C VOLUME. ANG CONTACT ANGLE AND TEST FOR 

?DS C CONVERGENCE 

207 2200 CONT TNUE 

?0B CALL STA 7E 1 1 Sft VF .USA VE * ZNE M.RS.ZS .RPS. ZPS* VS) 

209 VU = VUlL tvS.2S.2W AX ) * H5D .0 COO/ I V OL 

2S0 PHIF r D3L E I A IAN? I SNGL1 RPS ) .SNGl < ZP5 )) ) 

218 PHTC = OAT AN < R C ANP < 2 S I ) 

2)2 P N S T r PHIF - PHIC 

21 1 PH IT 0 = P HIT / ANUM 

219 3 F ff PHIF *1 T .0 .0000 .AND, P H I C • 6T .0 .0000 ) PHIl - PHI I « Z.ODnOePI 

285" PS f = C « V U~ UL P CT S / I 00 .0 0001 * •2 * 

2)6 * < 4PHITD-PHI3 )/ *30^0000 l*#2 

28 7 PS ST - PS I 

2)3 $COFO : ASAVE 

28 9 SGN s SS GN 

220 _ ^ WRITEIN01.2005 ) ft SA VE *P SI *ZNf W. VU .PH I TO 

228" If C PS I * GT . “ t PS C ) GO T 0 2000 ' ~ ' " 

222 C 

22 3 C 
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224 

225 
? ?5 
2 2 7 
72 3 
229 
233 
2 3 1 

232 

233 
? 3$ 
235 
235 
237 
233 

23 3 

24 3 
24 1 
2ft? 
2*1 3 
2 <1 VI 

245 
2*5 
24 7 
24 9 

24 9 
25D 

25 1 

252 

25 3 

254 

255 
255 
257 
259 
2 5 9 

253 
2fH 
262 

263 

264 

265 
2 5 5 

26 7 
253 
269 
? 73 
271 
27? 

27 3 

274 

275 
275 
27 7 
273 
27 9 


C END OF HUN tOOP 

c 
c 
c 

C CONVERGENCE ESTABLISHED* TABULATE IR*2) COORDINATES 

C OF THE FREE SURFACE SHAPE ANC GO R t AO DATA FOR 

C NEXT CASE 

C 

c 

CALL PAGE HU 

WRI I F I NO T *2006 1 1 0 A S H . I r 1 ,5U) *9 ONONO* PH I D *UL PC T • 10 A SH *1 = **5U ) 

T =0.0 DUO. 

DO 1100 I = 1 * J S A VE 

S OLNC I* 2* IS AV FT = ZNFW * S 0 L N \ I* 2 * TS AV FI 
CALL PRIKTS1I*NE0*ISA VET 
T = OFt T AT* DHL E (FLO AT (111 
1 l 00 CONTINUE 

GO TO 100!) 

C SURVEY COMPLETED* TABULATE RESULTS 

1 5 UO CONTINUE 

Hr? 

15 10 C0N7 TNU F 
12 r 11*4 

IfU? » GT D N A) 12 = N A 

CALL PA St HO 

WRIT F 5 NOT » 2H04I ( 0 AS H« 1=1* 1 001 « <SUPV PL If » I* M *:j=H * 12 y 

WRI ?E 1 NO T *20 Q 3 ) 
no a 5 2ti I = 2* N2 

WRI TE( NOT *20U9 1SUPVPL ( ? *1 *M*t SUR VP II T » J * 1 >* SURVPI 1 1 «J *2 )*J= I 1*12 J 
15 20 CONTINUE 
If r 12*1 

IF (If .GT. NA) GO TO 1 S4fJ 
GO TO 151 0 
15*0 CO NT TNU F 
C PLOT RESULTS 

CALL LOmOPL I5URVPL*N A- 1 .NZ-1 * SNGf <7 MAX). 51 ) 

30 TO IDQU 
C 

C ERROR EXITS 

C 

9000 CONTINUE 

WR TT f C NOT. 201 111 
STOP ' 

93 10 CONTINUE 

WR I IE I NO T .2020 ) 

ST OP 

3 D? 0 CONTINUE 

WRIT F( NOT. 2 0 30) 

STOP 
C 

C FORMAT SI# IF MF NTS 

C 

2 00 2 FORM# Tf^/4UX*?0.HRE SUL TS FOR X SWEEP 

1 IfVrM HX- INT ERCFPT *I2X* 5 HER ROR 

2 7X*I0HULLA6E PCI * 4X«I3H CONTACT ansif 

3 6X.MHFLIU0 ANGLE . 3X. 9HCAN ANRE 


• tt 
» 

* 

• / / 




of 





?3T) 
?9I 
?8? 
?83 
?9H 
?9*5 
235 
23? 
?B5 
?3 9 
?9Q 
?9l 
29? 
?91 
234 
295 
295 
297 
?99 

299 

300 
3CH 
30? 
in? 
10 6 ? 
505 
305 
307 
309 
309 
3 * 0 
3S 1 
3 1? 
3« 3 
3 1 *5 
31* 


NO SOtOf ION * *•** 


?srji FORM AT JSX'o Of 2o6» 53 X*25 «•*•** 

?003 FORMA H // ®*»0X oFOA L 

1 58X* ISHSURVEY SUMMARY 

2 30 X *5 mi 

3 I ?X o t NX * I OX * 5 < 4 H A “ *f6o?*9XI 
2003 f ORMA H l 5 X ®5HCOORO*5 <5 X *5 M VUPCT*?Xe6H 
2005 FORM AT S///2X* 3 1HMN IMUM ERROR SOLUTION 


| 22H TRUE 

I 1 2H T RU F 

l 2 2H TRUE 

? 2 2H TRUE 

* 7 

2306 format j i hi «/ /y • i nx «*n ai 


;RROR - 

X- INT €R CfPT = 
ULLAGE PCI - 
CO NT ACT ANG = 


? 

? 

3 

4 

5 

6 
? 
9 
9 


PHI T3 > 
FOR «COF 
tD I ? • 6 * t 

# oi i .«•■/ 

pD I 2*6 */ 
*• Dl 2 ®6 


SH* PE FOR 


T 

* / 

* / 

9 f t f 

I 

9 / > 

q* Cl ? 


, 6 1 3 H TS/ 


P R7 ME ?G X * 


I0X.35HFRFE SURFACE EQUILIBRIUM 
25X « I 5HM0NG NUMHtR - *OI?-5 

25 X* I 5HC0N1 AC T ANGLE - *0!?«6 

25X , I 5 HULL AGE VOLUME - *0<?»G 

J OX. 501 I 
7X,5HARC 

4X»6HLENGT4»9X*IHR*I?X»IHX*9X*7HR 
?HX PR IMtiSX.5 HVST AR 

> 

2Q07 FORM AT f ? <5X » D$ FI? .311 

2 009 F0RMATH0Xs r l0.3*5<5X*F6.7»2X«P6.M» 

2009 FORM AT I / // 2UX • 50 41 • / 

1 2 0 X *344 SOLUTION TRAJECTORY F OR ACOF - 

2 20X.S0AI./// 

5 1 r 

2313 F ORM AT < 3 6 H HALVING LOOPt EXECUTION TERMINATED 


p 3 t 2 a 6 / 


I 


•> / 

9 t 
9 t 
9 t 
! / 

9 t t f 
9 f 

9ft 


SOtUTJON FOR ACOF :»DI2o&9 


2011 FORMA Tl // /?X «35HN0 MINIMUM ERROR 
I J $H PLUS OR MINUS *DI2.5) 

?{I20 FORM# II S2H RM»X AND/OR XM*X NOT OF F INF 0 .EXE C lH I ON 

2D3Q FORMAT <51 H N6 OR NX GT 50 FOR S E AR CH -FjALa t» FMrUTION T ER W IN AT EOl 


END 


3PRf CAN 
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F M3 UGG*IN?07*TP FS. CAN 

I 30U9LE 


PRECISION FUNCTION C4N1Z) 


7 


IMPLICIT fSOUSLF PRECISION fA-H*0-Z> 


3 

4 

5 
S 
? 

3 

3 
I 0 
! 1 
* ? 

I 3 
i 4 
IF 
I S 
I 7 
1 3 
I * 
20 
?i 
2? 

23 

24 
75 

25 
2? 
23 
23 

?n 

31 

3? 

33 

33 

35 

36 

37 
3? 
3* 

40 

41 

4 ? 

43 

44 

45 

46 
4 7 
4 3 

4 ? 
50 
5t 
52 

5 3 

54 

55 



c EVALUATE. FOR A GIVEN Z. R OF TKf CONTAINER OR THE DERIVATIVE 

C OF R WITH RESPECT TO Z 0? THF t ON T A I NFiR 

ccccccccccccccc ecccccccccccccccccccccccccccccccccccccccr cccrcctcccccccc c 

cccccccccccccccecccccccccceccccccceccccccccccccccccccctrcccccccccccrcccr 


COMMON /CONSTS/ PI.ANUM 

COMMON /PAR AMS / ACOFt FCOF*.?M AX.RMAX 


FNTRV RCAN4Z) 
IZ = 1 

GO TO IOU 
ENTRY RC^NPfZ) 
IZ- ? 


c 

C CO OF. COMMON TO T HF FV ALU AT ION OF POT H R it I 
C 

1 31) CONTINUE 

I c l 7 .GT. *7 .3000) 30 TO MO 


JZ - 1 

X z 27.3)00/47.3300 
GO T C 199 

1 iU I M 7 . 6 T » 55 .9000) 50 TO »2U 
JZ - ? 

X : 3.3000/3 ,8)00 

GO TO 199 

120 I F ( 7 . G T ■ &S. 5)001 30 TO » 30 
JZ = 3 

X : * . 1 3 0 0 n . 5 1 0 0 

GO TO 199 

130 I F t ? .61. 123.2DOO) GO TO 140 

JZ = n: 

X : 32.2900 

GO TO 199 

140 I F 4 7 . 6 T.« 148.5900) 30 TO 400 

JZ r 5 

X r 13 2.2 ODD * * 2 ) - 132, 2300 / 2 3 . 30 00) * * ?♦ 1 Z- 1 2 3 . ?D00 1 * *2 
I FIX . GE • 0 ,0 000) X = OSQRTTX) 

ifcx .Li. n.uono) x = u.cmoa 

m CONUNUF 

SO TO 120 U. 300) • IZ 


c 

C CODE FOR Rt Z ) 

c 

200 CONTINUE 

2tQ* ?Z0. ? 30* Z41 • 2 501 » J7 
210 cV*F?Wx*Z 
return 

220 CAN - 27.3D0U ♦ X * I 2 -4 7 . 3 0 00 ) 
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56 


R FTURN 

57 

?30 

caw - 3 B o 1 DO 0 4 X* 1 Z -56 .9000 > 

5* 


RETURN 

59 

290 

caw r X 

60 


RFTURN 

51 

250 

caw r X 

6? 


RFTURN 

S3 

C 


69 

c 

CO OF FOR R PR TMF <Z) 

S 5 

c 


66 

300 

CONTfNUF 

67 


GO 101 3IO*3?tl. 330o340.350) • JZ 

63 

310 

CAN r X 

3? 


RFTURN 

70 

52 0 

CAN - X 

?1 


RETURN 

7? 

330 

CAN : X 

73 


RE T URN 

79 

340 

CAN - 13.0 030 

75 


RE 7 URN 

76 

350 

CAN : 14 73.697 7 000 - 3 .3 1 9 706 000 *?» / O. .0 00!) * tX ♦* .0 0- * « M 

77 


RETURN 

73 

C 


79 

c 

error fxits 

BO 

c 


9 1 

4 00 

CONI 7 NUE 

9? 


CAN 2 1 .00* It) 

93 


RE T URN 

99 

c 


95 

c 


36 


f NO 


9PRT SOrffiCK 
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FM3UG6HlN2t)7*TPf I* GO HACK 


) 


105 ICftL 

FUNCTION 50 MCKISOOD) 


> 


IMPLI OIT 

OOURLE PRECISION 1 


5 

* 

c 

l 

cccacccccccocccccccccctcccccctcccccccccctcccccccccccn crctccr cccccccccca 

t* 

3 

6 

y 

u 

C 

ROUT INF 

TO determine when to exit integration LOOP 

f 

3 

t 

c 

ROUT INf 

TO oetfrmine when to EXIT INTEGRATION 1 OOP 

3 

1 D 

L 

ccc acccccrccfccccccccccccccccccccrcc rcccc cccccccccccctt ccocccc tctcccccccc 

1 ( 

c 




1 2 


LOG! CAL 

GOOD 


1 ? 


COMMON 

/VECTOR/ Y*F)*V0TI5> 


1 3 


COMMON 

/ T IMESS/ OET T AT* T 


15 


COMMON 

/PARAMS/ A C 0 F *3C0F » ?MA X * RHA X 


1 5 

c 




77 


50 AACK 

: .TRUE. 


1 3 


GOOH - . 

TRJ f . 


13 


I F f VI 1 t 

.61* PMA X ) 50 3 A C K = .FALSE. 


211 


IFIY ( 1 ) 

.lt . n.unmii go fa cm = .false* 


2* 


IF 1 ¥ < 1 1 

.11. P. D3U0) E COD : .FALSE. 


2 ? 


I F(Y t7T 

.61 • n.ontin .and. vmi .lt • u.doddi 

GO RACK - .FALSE. 

2? 


IF< M 2 ) 

,GU 0*03 UCI .AND* VM) .Lt. tl.UOODl 

GOOD t .FAfSf. 

2M 


I F (V 1 

.LT . 0.UDU0 • AN 0 • V(3) * GT . 0 .0 OOU ) 

GO PACK - .FAtSE. 

25 


IF C VI 2 5 

.LT. O.tmin .AND. YI3> .5 T . 0.0000) 

G 0 03 - .FMSt . 

25 

c 




27 


RETURN 



23 


END 




*PRT L 0 *0 a L 
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FM3UGGHIN2D? *TP FS 
1 

2 C 

3 C 

ft C 

5 C 

C 

c 
c 
c 
c 


& 

7 
9 
3 

« n 

i ! 

8 i 
t? 

I as 
8 5 

8 G 
8 7 
I 9 
8 9 
?n 


, LO«OPL 
SUBROUTINE 7 


L090PI TO A TA *NA »NZ . ZMft X *KOft Tft > 


: INPUT 0670 ft RR 6 V 
= NO OF A CO F VALUES IN flAFA 
r NO OF 7 SWF?p VALUES IN DS 1A 
= MAX AXIAL TANK COOROINATF 
K0QT6 r ROW DIMENSION OP Oft 7 A ARRAY IN Cft Li I 

NOT F MAIN PROGRAM MUST CALL 10£NTt1> 

a N 5 :«LL ENDJ09 70 TERMINATE 


N ft 
NZ 
2 MAX 


NS PROG PA I* 

TO initialize 

SC 9 0 ? n 


SCftOZU 


DIMENSION 06 1ft I K D 6 T A .KOft 1 A *8 > e 1 Y1 82 ) e T X i I 2 &»lSYt50)»Y«S0 I 9 XI SO > 


HAT A 
0 ? H 6 


WOT 
1SY 
24 2 


t 

/ 

* ? H 


3*24 ft B ZH 5 * 24 6 *Z4 7.24 9*ZH 3 *?4 10 » 24 It *ZH 1 2 • ?H I 3 „ 
*248*3.2 48 5*2 Hi 6 * 2 M 7 s 2 Hi 8.? HI 3* 2 420 * 2 H21 *2 42 2*2 423 e 2 H? 4i* 2 425*2426* 

*2H27.2H28*?423®24 3D.243i»?45?»2H33*243 i «*2H35*PH35.2M3?*2H38*?4 339 
*24ftD* ? 4»l J * 2 Hft 2 * 2 Hft 3*2 Hft ft * 2 Hft5* ? HftG. 2 HU 7#2 Hft 3*2 4ft3* 2 HSi! / 

00 ? 1 - 1*12 


27 

2 ? 

23 

2 4 

25 

26 

27 

28 
23 
30 

3 1 
3? 
3 3 
3<3 
55 
36 
5? 
38 
33 
40 
■8 $ 
tg ? 

IS 3 
ft ft 


TXII) r SH 
2 T V i T> - 6 H 
C 

7 X <5 I : fiHl ANK A 
1 X< 6 ) r & H XI S IN 
1 X i 7) Z 6 HT FRCFP 
1 XI 9 > : SH T 
C 

NR : N7 ♦ 1 

N C - N A + I 
C 

C ULL. AG F VOLUME PLOT L z 8 
C CONTACT ANGLE PLOT L = 2 
C 

00 1 00 L= I *? 

GOTO ( I U I 9 I n 2 I 0 L 
SO f YHft X " 8 DU • 

7 Y<5| - SHULL AGF 

TV(S) - SH VOL UM 
T Y C 7 ? r SHF P F'R 
1 V8 8 > r 54CE NT 
GO 7 0 8 03 
&02 VMft X r 8 9 0 □ 

T Yd 5? z 64CCNT AC 


ft 5 7 VI S > r &H T 8 N3L 

46 1 Y 4 71 = SHF OE G ' 

4 7 1 VI 3 ? r SHRcES 

*8 103 CONTTNUt 

S3 3 £ 

5D SFR = 0 

SI HR I T E ( NO T . I DUO 1 I TYI II«lrS*9 >* f T XI 1) 9 f z 5 »S I 

5? 1 ODD FOR^ATU F 6 0 I0X.16HPL07 SYMBOLS FOR»Y« 5 X* A AS* ft H V S »ft^S.//* 

53 " a SOX. S4SYM40L® 1 5 X * ftHfiCOF*//) 

Sft c 

55 00 120 KA z ?*NC 



75 


56 

57 

59 
5? 

60 

5 t 

6 2 

5 3 
69 
55 
66 
57 

69 
53 

70 
77 
7 ? 
73 
79 

75 

76 

77 
73 
73 

90 

91 
8? 
93 
39 
95 

36 

37 
99 
93 
90 
31 
97 
93 
99 

95 

96 

97 
99 
93 

100 
10 ! 
ID? 
103 
109 
70S 
106 
t 0 7 
103 

6 09 

I I n 

I I I 


AC “ T) AT M I • K A* S ) 

IS - ISYTKA-I ! 

writ f< not * < non IS * AC 

!0U1 FORMA 1< I 2 X *4 Z * ! 0X*E I 7 0 B l 
C 

KNT = 0 

DO 30 KZ = ? * NR 
XX = 04 IHKZ.i r I ! 

I f ( L .FO. I ! YV - OATATKZ. KA.M 
IML „EQ. 2) Y Y - 36 74 I KZ *K4 *? J 
IFIYY .GF. 0.0 • AN G • YY .L F • YMAX) GO TO 31 

IF « KNT „E0, 0) 60 TO 3 0 
C 

75 IFITFR «EQo 0! CALL GUI K 3 L 1-1 *0 * * 7M4 X * 0 • • YK4 X *35 • T X * I Y * -KNT * X * Y ! 
IfdFR 0 FG. 19 CALL 0UIK3LC 0* 0 „* 7 N AX 1 0 .# Y M AX • 35* T X • T Y * - KNT • X « Y I 
C 

CALL XSCLV<<X< ! » * IX R AS* IX FRRI 
CALL YSCLVMYI I »*1 YRA S«I YE RR I 
CALL PR INTV <?* IS* IXR AS* IYR AS) 

CALL XSCL V! I X!KN 1 !.I XR4 S.I X£ RR) 

C0LL YSCLVt <V (KNT)* IYR AS* IV t'RRI 
CALL PRINTVl?#! S.I XRA 5*1 YRAS! 

C 

IF R = ! 

KNT =0 
SO TO 30 
31 KNT “ KNT ♦ t 
XTXWT) : XX 

Y C KNT ! - Y Y 

XFCKZ *E Ob NR'! EC TO 75 
30 CONTINUE 
120 CONTINUE 
100 CONTTNUF 
C 

C CROSS PLOTS ULLAGE. VOLUME VS CONTACT ANGLF ONF/FRAMF 

c 

I Y ( 5 ! = 6 HULL A GF 

TY157 r 6 H V0LU3 - 

TH?) : 6HF PER 0 

T Y I 9 > = SHCt NT 
TX <51 = 6HC0NT AC 

1X16) = 6H 7 4 N3 L 
T.xm r 6HF Of G 
1X13! = 6 H REES 
C 

WRI TO NOT *>00? ) I TV! I ! *1 = 5 *3 ) • < 1 X 1 1 ) • 1= 5 .9 I 
< 302 FQRMAIftHl* !0X» I6HPL0T SYMBOLS FOP*/* 5X*3A6* 9 H V S * A A6. / /* 

* BOX* SBSYMBOL* »5X* 9HJNTE RCEPT* //) 

C 

DO ! 30 K Z= 2 *NR 
IS = ISYIKZ-U 
ZI = 0 A TM KZ * 1 •! ! 


1 50 U fi ll F < NOT • 1 U0 1 I IS * Z I 




76 


a $ ? 

1 I 5 
9 i a 

9 1 5 
a 9 & 

0 7 
t \ 3 
6 59 
4 20 
52 4 

1 2 ? 
12? 
124 
I 25 

4 2G 

5 2? 

1 29 
523 
9 30 
9 34 
I 3 2 
B 33 
9 33 
435 
I 3S 
5 3? 

9 39 
5 33 
1 9 D 


K NT =0 

DO ISO KZ: ?*NR 
X& = na? At KZ* K Ae ? 5 
7V i DA 1A J KZ at A • > > 

If<VV «GF. 0,0 o 0NCj* VV -L4F. 100,01 00 10 1*51 
IF 6 KNT «£ Q, 0> GO TO 450 

675 JF« SFR ,E 0 o D> CALL 0U1 K 3 1 1 - B * 0 . • » 9 0. »E1. ♦ I 00 . t35 *1 X • T » *-K N1 • X • Y J 
IFC1FR oFQo 4! CALL GUIK3L< Q«tl.fl9n»*D.*ir4n.B35«TX»TT»“KNT*X*YI 

00 i GO Kt=4 ®KNT 

IS - I $Y 4 3SL < Kt ~ * > 

CALL XSCL’VMX IKt B • IXR AS* IX ERR! 

CALL YSCl Vl I Y4KL B*I YRA S*l VERR B 
9 GO CALL PR INT V ( 2 « IS« TX R AS • IYR AS 5 

IFR - 5 
KNT - 0 
GOTO 950 
15 9 KNT " KN T * 9 

I F ( K ftl? 0 FQ o 9 4 ISL = KZ - 4 
MINU : XX 

V4KN75 - 77 

IF 4 K 7 o E Q b N R 9 30 TO 175 
J 50 COM INUf 
99D CONTINUE 

RE T URN 
END 


3PR? PASEHO 
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F M3 UG GRIN 2D 7* TP FI, P AGf MO 


I 

7 

3 

ft 

5 

S 

7 

3 

3 

i n 
1 1 
« ? 
13 
I 4 
IS 
I S 
1 7 
I 3 
1 3 
2D 
21 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 


SU3ROUUfcfc P ft 6 E HO 

COMMON/ LSI &RT/ IRUNNO* 10 61 F» NP 6GF* UN AME< 3) « T IT I El (I ?!* 1 HI E2 M ? | 
DATA N I ■? *N0 T /5 * S / 

3 RINGS UP NEN PAGE ftNO PUIS HEADING IT TOP. 

INTERNAL VARIABLES. I IRftNSF F RRE-0 THRU COMMON 1 !. 

IRU NNO = RUN NUM HFR * I AG FORMAT) 

IDA TE r DATE, (AS F ORM A T ) 

NPftGE. r P*GF NUMBER. 

UNAME r USE*RS NAME. 13A6 FORMAT) 

TITLE! - FIRST TITLE. ( I 2 66 FORMA?) 

IIUE2 - SECOND TITLE. II2AS e 0 RM A 1 ) 

2001 FORMAT t 3 H ) R UN NO. *«6*42X*5X .6 X *M ? X «SH PA GE NO . .1 E / 

♦ SSX.R HR UN WY* 1 X * 3 66/ I OX* 1 2 66/IUX. * 2 66J 

NP AGE - NP A G£ ♦ I 

NR I TEINOT *2DU t > I PUN NO *NP A $E • UNA ME * 1 1 U E ) * TT TL E 2 

RETURN 

EN3 


3PRT PRINTS 




78 


FH 8 UGS«TN 207 

I 

? 

1 

«4 

*: 

G 

7 

9 

=5 

10 
I 9 
> ? 

8 J 
I * 


* TPF S*PRS NTS 

SUBROUTINE P« {NTS CN.N FQo I) 

IMPLICIT DOUBLE PRECISION H-H*0-?> 



cc erect cccccccccccccceccccccccccceccccccccccccccccccccucccrxccccccccccrt 


c OUTPUT ROUTINE - PRINT SUTf c OR EACH fULUT 



c 

COMMON /TIMtSS/ D£lHT»l 

COMMON SOLNtlOUQ'S'Z) 


C 


a s 


rm a NOT /£/ 


IS c 


8 7 T S “ SNGUl I 

IS WRITE* NOT »IQUO i IS* « SOI M N *J *T I • J = I * NE Q > 

89 ^ 3 nn f ORM 6T < FID » 3«5 01 T ,S) 

?D C 

2 8 RETURN 

2? I N3 


5PRT RCSN2 



79 


FHBUSGHIN207 * 1PF S*RCA N? 

| DOUBLE PR EC IS ION FUNCT ION RCAN?f2» 1 

? IMPLICIT DOUBLE PRECISION ia-h.o-zi { 

5 C ' 

% ccccccccc cccccccccccccccceccccc tcccccccc cccc cccc cc c-ccetccctcccccccccccc ci 

sc 

5 C ROUTINE TO EVALUATE 1HF SQUARF O p R OF THE CONTAINER TO *E l 

Y C use'll IN COMPUT IN6 THE ULLAGE VOLUME 

3 C ' J 

3 ccc cccccccr ccccccccc cccc cccccccotccccc cccc cccccccceaccccccccccc ccccccccc 

10 c 

i I RCAN? = CRCANIZI I ♦*? 

» ? c 1 

13 RETURN 

TH END 1 


4PRT RUNK1A 
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FH3UGG«IN207«>?PFS o RUNKT « 


SUBROUTINE RUMM IS ? N t OTNoMI) 

IMPLICIT UOUHLf PRECISION <A~H«G-ZI 


5 

9 

c 

cccccccccccccccccccccccccccccccccccccccccccccccccccocccccfcrcccccccccccr 

5 

s 

L 

c 

RUNGF-KUTT a- GILL NUMERICAL INTEGRATION ALGORITHM 

1 

7 

9 

c 

cccoccccccccccccccccccccccccccccccccrcccccctcccccccciccoccccf cccccccrcccc 

9 

s n 


CO WH ON 

4 QPRKT A/ QR K IS 1 

« PR K 1 9 J 

» } 


CO*§HON 

/TI*E SS/ DEL U T 

9 T 

i ? 


COMMON 

/ V FCTOR/ Y ) * V OT (5S 

1 3 

c 




j $, 


00 1 20 J 

- I 0 9 


IS 


JTL r J 



i s 


DO M □ 

I z IoNEOTN 


17 


? - YO 1(1 } OEL TS T 


1 8 


GO T 0 

I 111 3« 1U I • 19 i » 1 US) * 

JU 

53 

6 OS 

R r 

PRK S JJ L 1 *< 7 - QR* 1 1 

V 1 

20 


GO 

TO 10 ? 


21 

1 03 

R r 

PRK I J1L 1*2 - OR* II > 


2 ? 


GO 

TO 10 7 


23 

1 OS 

R - 

tl - 2.0300*GR< * I 1 » 

/ 6.00 00 

29 

10 7 

Y < IS 

: YU1 ♦ R 


25 

1 ID 

0 RK II ) 

z ORK U ) * 3 rOOOO*R 

- PRK t J1 L )* Z 

2 S 


IF < JIL. 

6 Q o « e 0R. JK.tO.31 

T z T ♦ OFLT AT/2 .0 D30 

27 

S 20 

CALL YO C I 



? 9 

c 




29 


NT r NT * 

) 


3D 


T = flPLflFLOAT <NT))*0£ll ST 


31 

c 




3? 


RETURN 



33 


EM3 




3PRT SAVE 







FMBUGG»lN?n7*TPFS*Sfc VF 

i SUflROUY INF S AV F T N R 0 Id * NFQ* I A J 

7 IMPLICIT D OUT L E PRECISION M-H*0-Z) 

3 C 

* ccceecccccrcccccccccccccccccrcccccceccccctcccccccccrccrccccrceiccccccccc 

s c 

5 C ROU TIKE TO Sft Vf INTEGRA TEH STATE SPACE SOLUTION 

T C 

5 CCCCCCCCCCCCCCCCCCCeCrCCCCtCCCCCCCCrcCCCCCCCCCCCCCCCCCtCtCf-LLCtCCCCCCCCr 

c 

ID COMMON SOIN I » UOU*5 *7 I 

II COMMON /VECTOR/ YlSltVCTlSI 

»? c 

13 00 10 l:.I.NFQ 

m 1 C» SQL M N RO W * I • 1 A I r yin 

IS c 

1= RETURN 

» 7 FNC 


*PRT START 
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FM8UGG RIN2D7 *TPF»- START 

1 SURROUT INF ST ART 

? DIMENSION MON TMN I I 2 > *MONTML I T ? I tA 3* PA II 2 l > 

} COMMON /L5T ART / TRUNNO • 10 AT E* NP 4GF*UN AMt <31 « 1 Mi F» « 1 7 I • T 1 1 1 E 2 f» 2 > 

4 DATA Nil * NOT /5 *o / 

5 DATA MONThN/?HfU*?Hn?#2HO^*2MO ? l#2HOSt2MOo* 

5* * ?4U7 t? s 4{i3»2rl03 ■2H(Q«2Hl »«2Ht2/ • 


7 
B 
3 
It) 
i i 
i ? 
i i 
1 1 
i s 
IS 
I 7 
» B 
I 3 
20 
21 
? ? 
23 
2* 
25 
? 5 
27 


* MONT HL t 2 HJ A» 2 HFfi .2 HMR* 2 H AP • 2 HMY • 2 HUN* 

* 2HJL flHtk U*2HSF »?HCC *2HN0* 2HDE / 

OAT A It ST t 0 ✓ 

c 

13 01 FORM AT C A6. BY T AG I 
1 00 ? FORMA T * I 2AS ) 

2Q03 FORMAT (3SHIFNL OF INPUT OAT A NAS HkfeN R t ACHE 0* ) 
C 

IF (f151.fQ»0t CALL IOfeNT < 3# A OAR AT I 
1 1 S 1 " I 

READ (NIT'IOUM IPUNNOtUNAME 
IF IIRUNNO • NE . BHSTCF) SO TO 10 
CALL FN 0 JO ** 

WRITE f NO It? HU 3) 

STOP 

C 

13 READ < NIT .1002) TITLE! 

READ INIT.JDOZT TITLE2 
NP AF F - 0 
RETURN 
‘‘i FNC 


3PRT ST ATE 

o 





FM3UG6MTN207*TPF$*>ST AT F 


1 

I 0 
\ 1 
I 2 
I 3 


SU* ROUTINE S TA TF IT * J»ZT *R S *ZS «-RPS*ZPS • VS > 
IMPLICIT OGUHLF PRECISION < A- H* 0- 7 > 


C 

c 


HAVING 8R AC Kf T F 0 CAN# FINO A GOO 0 ST AT F 


ccc occccccccc CCCf ccc ccccccccccccc CCCCCCCCCCCCCCCCCCGtCCCCCCCCCCCCCCCCCCC 


COMMON S OL N < 1 flUDt 5* l 1 

CALC* K 0 7) I SOLNIJ- I . K* II ♦< FACT OR* IS Of. N«j* Kf IT- SOL NfJ- «# K« It I > *Z 


i 4 
15 
i fi 
17 

1 1 
1 1 
20 
21 

2 2 
23 
23 

25 

26 
27 
21 
21 
ID 


R C IN - RCAN(SOLN(J- 2f II *Z I) 
RCOUT i RCANI 50LNI J* 2*1 >♦ 71 > 
Ft - RCTN - SOi.NI J- 1. 1 « 1) 

F 2 r SOL N ? J • t * I > - RCOUT 

F ACT OP = Ft / I FI ♦F2I 
RS ' CALCIttO.OIDOl 
ZS “ C AL 0( 2 • z n 
RS = C ALC t » *U.UDUn ) 

ZS - CALCIZ.ZIt 

rps - calc *3 .n.cionin 
Zps - caloi I* o .0 orn t 

VS t CALC* SfO.aDflO) 

RS$ r OS GRT f RPS * * ? ♦ ZPS**?> 

RPS : RPS / RS S 

ZPS - ZPS / RSS 

RETURN 

FNO . 


}'i 


^PRT SIMPS 
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FM3UGGHIN207*?PFS. SIMPS 

1 D0U3LF PRECISION FUNCTION SIMPSIA.3.F) 

5 IMPLICIT OOUMLF PRECISION <A-HiO-Z) 

^ cccccccccccccccccccccccccccccccccccccccccccocccccccixcceccccccccccccccci 


S 

7 

3 

9 

i n 
it 
i 2 
»3 
J * 
15 
1 s 

1 7 
« 9 
13 

20 
? s 
22 
23 
23 

25 

26 
2? 

29 
23 

30 
3 » 
3 ? 
33 
39 

35 

36 

37 
39 
33 
9 0 
9 ! 
9 2 


C FUNCTION TO NUMERICALLY INT EGR AT E THE FUNCTION F 

C FROM A TO 8 USING SIMPSONS RUlF 

cccctcccccccccccccccccecccccccccccccccctccctcccccccccccccccccccccccccccc 

c 

DMA I N T / 9 D i 

c 

c A-LOUFR LIMIT OF INTEGRATION 

c fl-UPPER LIMIT OF INTEGRATION 

c F-IMTFGRANO FUNCTION IDEC L ARID EXTERNAL IN C A U I N6 PGM) 

c IP AR- P flR AMFT tR PASSED TO INTEGRAND FUNCTION 

C 

c INITIALIZE PARAMETERS 

C 

TNOH = ( M- AT / net M FLO AT 4 INI M 
H = T WOH / 2 ,03 DO 
SUNFNnrO .0 uni) 

summio = o.oano 

c 

C TWOH-TNTERVAL 

C H- HAL F INTFRV AL 

C SUMMENO-SUM O r FIX SUB I)* FOR EVEN I 

C SUMM10-SUW OF FIX SUB lit FOR ODD I. 

C 

C EVALUATE S SJ ME N 0 AN 0 SUMMIO. 

C 

DO 4 K z i . INT 

XiA+DBLE f e LOATIK-T )) *1W0B 
SUMENO = SUMENO * F TX T 
I SU9MID = SUMMIO ♦ ? I X*M ) 

C 

C RETURN ESTIMATED VALUE OF THE I N TF?G RA ND 

C 

SIMP5 ” I 2.H0D0* SUME ND*3. 0000 * SUHMI D-F I A J+FIB » ) *H/3.O0UO 

C 

RETURN 

END 


3PRT VULL 
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FMBUGGBIN?07*TPF S.VULt 


I 0 0 U*t f PR E CIS ION RUN CT I ON VUt L I V St *3 A* 2 W AX 1 

? IMPLICIT DiOUSLE PRECISION U-HfO-Z)__ 

"5 c " "" : " - 

A cccccccccccecccccccccccccccccccccccccccccccc cccccccccccccj^cccccccjcccccce 

5 ~ c ' : ' " 

3 C COMPUTE ULLAGE VOLUME FOR A GIVEN SOLUTION 

7 ' C 

s cccccccccccccccccccccccccccccccccccccccccccccccccccccccccecccecccccj;c«c 

3 C 

10 _ EXTERNAL __ RC A N? __ . _ 

I i " ’ COMMON / CONSTS / PI* ANU M 

I? c 

13 VULL - P T* < VST AR ♦ S IMPS 4Z A.ZWAX #R CAN?1 I 

M ... C ..... ..... . 

1 5 RETURN 

I S E ND 


iPRT tOOT 
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F MB UGG8IN237 * tPFS, V OOT 


I SUBROUTINE YDOT " - -- 1 

? IMPLICIT HO U 8L E PRECISION 

? c • " 

t» CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCGCCCCCCCCCCCCCCCCCCCC 

~ 5” C “ * : i 

6 C ROUUNE TO COMPUTE VOOT « A FUNCTION OF V 

? c " ' 

8 cccccccccccecccccccccccccccccccccccccccccc cecccccccctcccccccccccccccccccc 

9 C i 

in COMMON /VECTOR/ YISI.YOTISI 

- ri ~ COMMON 7P V RAM S / AC’OF tffCOF t ZMI XTRH* X *“ I 

1 2 C 

I? R55 t 0SQRTIVI3)Vyi3) 4 Y|tn#V|«IM 

l«* V C 3» = vm / RSS 

15 tm : y € ^ > 7 RSS 

$6 CCOF = D.DDOO 

I T "I F t Y1 I T . WIV 0 i 00 00 T C C OF "z V ITT 7 ' V I I I "" ""7* 1 

19 TRM - ACOF ♦ 8COF*Y<2> - CCOF 

1 9 C 

2D YDTIII = V f 31 

21 VOT 121 : fji I ” " i 

2 2 VDTU1 - - V TB > ♦ T RM __ • _ 

23 ~ YD T I li ~z~ Y 15 M TRM ~ ' ~ - 

2 V OT < 5 1 z Y < 1 1 *Y < I I *V (4 I 

"25 ‘ C " 1 

2* RETURN 

2 7 t N O 


iF IN 



6.1.2 Listing - Vibr a tion Analysis C om puter P r ogram 
and Associated Subroutines 


3 H I L I P B I N 2 0 7 * F 2 » rt A I M 
?. C 

1 C— TTAIN ~P R"tTGTvA7T TO CTTl TrTnrXTF l TT7 — BONU N U M fi f ST T (TfTK S L CTSK TTGOE^ - ' 

M C DEVELOPED t > V K L A 0 H L F_ N , F E B K n A R Y i 9 7 5 • 


COMMON / RAP2 / I <(i H D 2 
T. (JTi H O'N / ;\A P5 T i A F< DT 


~C“ TTTP'U T D aT a' PE A l> "I M TR t S PROG h A M •' 

c jo call start 


FORMAT (ZA6J 

C Call DGEN> (SEE SUBROUTINE for INPUT) 

T - FioF T “ “ ~~ “ • fIj rTTaTt Vh7>) 

C If < H OPT , E G ® 6H M 0 D ED ) CALI M 0 D E 0 (MO INPUT R£ Q HIRED) 

T TV (MUPT • E~G . 6HU0UO } c7UT“Htf DES (SEE SUB^T FUR TTTFuTT 

c i o p t format ( a 6 j 


uhei 


M 2 

10 call start 


' t nr ~ 

M M 

Read ln i t» i ?xtd ~tf t~nt t » t a etd 

if (ifinit .eq. bh i n i t i l > call intape ( uhsvt i , tape id j 


MS C 

Mi 

rewind hotel 



D N U T X Y 

Call lbdgek* inutel .nutxyz > 

REWIND nuTel “ " ' — — — 

Call LBF I NE ( NUTEL »NUT*YZ t NUTM *NUTK , Nl.'TLT , NUTST , no tux ,hUTKX t NUT&X , 
* NIITB I »NUTB2 ,i.jUTB3) “ — — ' 

Read init,ioiu> mopt 


6HM00ED ) 

GO TO SO 



6HM00ES ) 

GO TO 7 £ 



N E R R 0 R a j 












3 H I L I PB 1 H20?*F l , FINEL8 

_ r~ r : O^P lXEl?^ OT~=n^, TO!T!T^^r^j T 

j2 SUBROUTINE FIMElB ( X Y t , JD HE , E UL , !\jU T E L , N J , 

~ 3 • * ‘ ?^u^tiTnutr~ ’ T^rr^cn 

** * K«X • KRj , KR£ ,'NUTMX , f.'UTKX ,NUTi, NU T2, NuT3 > 


» UMfcNiiO'M xr4<KHA,|)| JpOFIKRJ»l># E U L ( K R E , 5 ) , VCU, Lv(J) 

i DIMENSION 71 ( 2R ,2'U , 2 ( 2 T , 2H ) j *3 12*1,2*1) 

DATA Kit '7 2 H / , 1 B LAN K / 6K ” / , I 1 / 1 7” ’ ~~ 

data n i r , m o t / 5 , 6 / 




SUBROUTINE TO CALCULATE (ON OPTION) FINITE ELEMENT,.. 


NibtHBLtU MASS HATKIa (ON UUTMJ, 

ASSE«OLE.O STIFFNESS MaTrIx (ON N U T K ) » 

"TV E C GIVES ELE FTE'T'TT D~UT TN TO 0 L 0 0 A L OOF . c XA NPLE 5 . . . " " 

lVEC(6)~H3H PLACES ELEMENT DOF 6 INTO GLOBAL OOF 33'** 

"' I vec 1 3T»H uni ts elemen t" oof“3 from Global "ooFr Tins cunstra r ns 

ELEMENT G 0 F 3 TO 7 E R 0 MOTION* 


DA I A ARRANGEMENT 0 M MUTM, m U T K F 0 R TH E ASSEMBLED MATRICES IS I 
SPARSE ( Y ) FOP h a SUBROUTINE FORMAT* 

'DA I A ARRAN G E’MENT CTiT” N U T L T , N U T K X , N U f S T , M U T B FOR EACH FI 7HTT“ 
ELEMENT { W ft i T T E N I fi SUBROUTINE FLUID, t T C I IS < « a O 
tfRITE (NUT ft ) N A H E *”* N EL , N H~~, N C , N A M £ L , I I B L A N K ,T= t , S ) , 

( ( vv { I f J ) , J = J ,.VR ) , J = 1 ,'C ) , (I VEC (I)*I-) t NC) 


N A M E L = FLU ID, ETC. 

LAST RECORD (TO DENOTE TERMINATION) IS, 

Wkite (nutw) I3l7an"k7Ti i , i«j , 33 > 

the following utility tapes oSE basic fortran reao, write* oo not 

use these tapes "In spas s e ’ iT\ f o r m a s u b r o u tines ,v n i c h u s e’fWil a 

SUBROUTINES yin, Your < because they use buffer in# buffer OUT 1 * 


N U U'lA, N U T K X • 

THE FOLLOWING UTILITY TAPES USE FORMA YIN, Y 0 U T * 

" WfM » nDTK » NU T i , n U T 2 , N J T 3 , " ~ 

CALLS FORMa SUBROUTINES FLUID , G« A V T Y , P A g£hO , SU R f TN , Y H V AO 2 , Y L E RQ i 

zzbqmb. 

DEVELOPED BY ™a BENFIELO, CS bOOLEY, RL 'AOMLEN. JANUARY 1973* 







[i] 


5 & 

c 



X , Y , Z LOCATIONS RESPECTIVELY. 5 J ZE (i\|j ,3 ) . 

s~r 

c 



HAY B F TQ U 1 V 'A 'L 1 5.XID T 0 V'TTJ IN CALLING PKUGRA~H. 

SB 

c 

JD0F 

33 

matrix OF JOINT GLOBAL degrees of freedom. ROWS CORRESPOND 


c 



'TO J 0 T ITT “N U M 0 £ R 9 . COLUMNS 1,2,3 CORRESPOND TO THE JOINT' 

60 

c 



TRANSLATION OOFS and COLUMNS 9 , e ,6 CORRESPOND TO THE JOINT 

6 J 

“ C 



- „v "-'‘tv' v-v “3- LV-w - . - „V-V,s rVi v- : ~ ' *#■' */■ 

62 

c 



HAY E£ EQUIVALENCES TO LV()> IN CALLING PROGRAM. 

63 

c 

tUL 

= 

HA IRIX Of JqTnT EULER Tn^LES ( DEGREES). RQ A S CORRtSPQND 

69 

c 



TO JOINT NUMBERS* COLUMNS 1,2,3 CORRESPOND TO THE 

6 S 

r 



global x , y » z permutation, s I z e { n j , 3 ) • may se 

66 

c 



E«(UI VALENCEO TO V (KRX* ( XY2' col r> I M ) ♦ n IN CALLING PROGRAM, 

“ — 57 

— r 

N LH"E"~ 

s 

LOGICAL NUMBER OF TAPE CONTAINING ELEMENT INPUT DATA FOR 

68 

c 



THIS SUBROUTINE AND' SUBROUTINES AXIAL. ETC GIVEN BY NAMEl. 


c 



IT NUT EL = 5, D A T A .<ILL BE READ FROM C ARDS. 

70 

c 

N J 

a 

NUMBER OF JOINTS OR ROWS IN MATRICES (XYZ), (JDOFJ, (EUL). 

n — ~ 

c 

— NUT J*", 

X 

LOoKaL H UMfitTR Of UTILITY TAPE ON WHICH ASSEMBLED 

72 

c 



MASS MATRIX IS OUTPUT IN SPARSE NOTATION. 

73 

c 



N U 1 H MAY be ZtK 0 IF Ha”S"s" MATRI X IS N 0 T ~ F G R M E 0~. 

79 

c 



USES FORMA YIN, Y 0 U T . 

_ 

TT 

T7UTK 

= 

logical Dumber of utility tape on which assembled 

76 

c 



STIFFNESS matrix IS OUTPUT IN SPARSE NOTATION. 

~n 

c 



r. U 1 K MAY BE ZERO IF STIFFNESS MATRIX IS N 0 I FORMED* 

7a 

c 



USES FORMA YIN, YOUT. 

79 

— r 

V 

= 

VECTOR CORK 5 p A C E , 

80 

c 

LV 

s 

VECTOR WORK SPACE. 

a I 




dimension size of v", lv in calling program. 

82 




ROW DIMENSION Q F X Y 2 IN CALLING PROGRAM. 




RoiT. DIMENSION UF JO OF IN CALLING PROGRAM. 

8'+ 




ROW 0 I M E i«! S I 0 N OF EUL IN CALLING PROGRAM, 

83 




lUuicaL nuMs'Er of util'ity tape on which element 

8 6 

C 



MASS MATRICES AND IVECS ARE STORED. 

8 7 

V- 



NUTMX MAY BE ZERO IF MASS MATRIX IS NOT FORMED. 

8 5 

c 



USES FORTRAN R E A D , .-.RITE. 

ST 

c 

N U T K X — 

a 

LOGICAL NUMBER OF UTILITY TAPE ON WHICH ELEMENT 

90 

c. 



STIFFNESS MATRICES (SAME AS GLOBAL LOADS TRANSFORMATION 

9 J 

c. 



MATRICES) AND IVECS ARE STORED, 

92 

c 



Nu'TkX may be zero if stiffness matrix is not formed. 

93 

— ~T 



USES FORTRAN READ, ft‘ R I T E • 

99 

c 

NUT 1 

= 

LOGICAL NUMBER OF UTILITY TaF-’E. USES FORMA YIN, YOUT. 

73 

r 

NUT. 2 "' 

S 

L 0 GiCAL number of UTILITY TAPE. USES forma YIN, YOUT. 

96. 

c 

N U T 3 

a 

LOGICAL NUMBER OF UTILITY TAPE. USES FORMA YIN, YOUT. 

O 1 

C f 

98 

c 

1001 FORMAT (66) 

~ 9"9 


2 0l' i FiJWMAl — { / / 9 ] a UTFTJ 0 I H T DATA U 5 1 (j IN SUBROUTINE FINED 

100 


2202 FORMAT {//3SX 97HJ0JNT data USED in SUBROUTINE F 1 NEL (CONTINUED)! 

■ roi ^ 


2U33 FqRHaT ( / f p X J 8 H DEGREES QF FREEDOM 

102 


* 


ibx ?8hglodal cartesian coordinates 

■nw 


* 


J Z x 2 2 H E 0 L E << A N G L E S fD EG REES > 

1 09 




/MX 1 I H T R A N S L A T T 0 N BX SHROTaTION 

~ 


tsr 


/ 2XSH JOINT 6X|HU 5 X 1 H V 5X“"1h^ 5X1 HP SX1HQ SXTrtR 

1 06 


. at 


UXIHX 11X1 HY i | X 1 H L 19X1 HX i £ X 1 H Y 1 2 X 1 H Z /I 

107 


2TT9 format (lx is, 3X 6i6, 3x 3F12.9, 9X ^fiu^) 

1 08 

c 




1 t) 9 “ 


r r 

t N U T M X * (1 T » (TJ R E n I N 0 :v U T M X 

l 10 


IE 

(NUTKX * G T 0 Cl REWIND N U T K X 

1 1 1 


~WfLT 

a -r 


ORIGINAL PAGE IS 
OF POOR QUAUTX! 












0 0 3 b j = i , b 

rp iJoor<i,Ji — . gT. n o ory -i j d u f = j l> o f < i 

35 CONTINUE 


PRINT JOINT OOF » Ail C 0 0 R 0 { N TPS, EULER ANGLES 


1 29 
1 30 

IF INLINE ,LE, 92) GO To <4 0 
CALL PAGE HD 


131 

132 

WrI TET l NOT , 20.32 ) 
i»R I. TE ( NOT , 2UG3 ) 


133 
1 3M 

NL I TIE = i 

93. WRITE (NOT, 2009) IJ, (JDOF(IJ,J), J=!*6), (XY2(IJ t J), J=l,3) 

i 


READ FINITE EI.E IWTYPT, 

so Read ( motel, i son nanei 

" IF (HAMEL .EL* 6 H RE T U ,R N ) q 0 TO 50 0" 
IF (NAMEL . E (3 * 6 H F L U I 0 ) GO TO IS l 


VTY) Q 0 TO 17! 
IF (NAMEL .E'W. 4HSUSFTM) GO TO IV? 


"NE R R 0 R ■ 1 

GO TO 999 


151 call fluid (xy2,jdqf ,eul.nutel,nj. 


* NUTMXiNUTKX, NUTLT,mUT5T 

* ill 1 , ft ? t w 3 , K R X , K 8 J , K R E , K i» ) 

'GO TO 50 

gravity element, 

r 7 I t A Cl G R A v T Y ( XTZTJ’iJoF . EUL,MUTtL,NJ, 

* N U T K X » 


* '«1»W2,W3,KRX*KRJ,KRE,KW) 
G 0 TO 53 

'SWTACt TENSION ELEMFMT. " ! ' ' ““ 

190 call surftn uyz , jogk ,eul ,nutel ,nj , 

__ nutkx , 

* «l l A2*W3 t KRX < KRJ t KRE,KW) 


0 5 0 


,NI DATA UN J I U K 

5 03 IF IMUTMX ♦ G T # 0) WrE (NUTMXI (BLANK* ( II I 1*1 *30-.) 


F INUTKX » G T * 0) ft RITE tfjUTKX) (BLANK, I I 1 » I s I ,30} 


A T R I c E S 
LL yz fr 
LL YTCK 












! 74 


999 CALL ZZBOMB < 6HF I N ELB . N£R ROR ) 




’’H I L I p 0 I >\Z 

37 * F 1 

* K 25 T i 



~T~ 



LUMP J L t K { X n s 1 } , ( E N U l V » C M N ) 



2 



SUBROUTINE K2ST1 ( X 7 , X 3 , Y 3 , S T , Z , K l ) 



3 

4 

C 


DIMENSION 1 [TT> l 5 - 


- - 

5 

t 

3 U b K 0 U r I N L 10 C A L COL A T £ F I N ITt EL E ME N T , , . 



6 

r 


STIFFNESS MATRIX, 



7 

c 

PUK a SUHtACE TENSION T .4 I A N OLE ELEMENT U J t H UnRES Tr A I NED BOUNDARIES. 

0 

r 

linear displacement melt is used* 



7“ 

C 

STlFTNbSS i^TATrnr^IS IN LOCAL CO'ORD INaTE' sYSTEFT 

» 


1 0 

C 

THE LOCAL COORDINATE SYSTEM ASSUMES THE PLATE 

TO LIE 

'IN AN X-Y PLANE 

i r 

C 

•UIM JUUMI 1 AT the x-y origin, joint ? Lies along the positive 

l 2 

c 

X 

AXIS, AND JOINT 3 IS IN THE POSITIVE Y DIRECTION. 


1 3 

c 

LOCAL coordinate order is 



1 4 

c 


0 l I * D l 2 , D 2 3 . 



T 5 

c 

ii HERE 02 IS TRANSLATION (OUT OF PAPER). 



1 6 

c 

DEVELOPED) BY R L H 0 H L E N • FEBRUARY 5975, 



1 7 

c 





1 3 

c 


SUBROUTINE arguments 



1 9 

c 

XT 

- I iM P U T LOCAL X COORDINATE LOCATION OF 

JOINT 

2 • 

2G 

c 

X3 

- INPUT LOCAL X COORDINATE LOCATION OF 

JO I NT 

3 * 

2 l 

c 

TT 

- input local y coordinate location of 

JOINT 

3 . 

22 

c 

ST 

a INPUT SURFACE TENSION ( FOR C £ / LF NG T H ) 

• 


23 

c 

L 

* OUTPUT STIFFNESS MATRIX. SIZE (3,3). 



24 

c 

<2 

a INPUT ROW 0 1 MENS I ON OF 2 IN CALLING 

PROGRAM. M 1 N» 3 » 

25 

c 





26 



A a X2*Y3/2, 



27 



CONST = ST/(4.*A) 



23 



X3MX2 * X 3 - X 2 



29 



Z ( I , l J s CONST * ( Y 3 * * 2 + ~T 3 MX 2 * * 2 ) 



3 2 



2 < 1 , 2 > --CONST * ( Y 3 * * 2 * X 3 * X 3 M X 2 ) 



3 i 



2(),3) = CONST * X 2 * X 3 H X ? 



32 



2(2,2) = CONST * ( Y 3 + * 2 + X3**2) 



3 3 



2 ( 2 , 3 ) = -C ONST * X 2 * X 3 



34 



2(3,3) * CONST * X2««2 



35 

c 





3 6 

c 

SYMMETRIZE LUtfER half. 



37 



DO IS J=l,3 



39 



Oo ID I *J ,3 



39 


TV 

2(3 , j) a Z ( J 1 1 ) 



40 

c 





4 i 



RETURN ~ ~ ' “ 



42 



END 




5PRT Fl.LBDGEN 




PHILIPBIN2O7*F».l-B0<3EN 


v. U i 1 > i Lu i\ l a H “ 1 l t 1 L Itf U l v 55 L H f 

subroutine l b d g e n inutel.nutxyzi 
t o M m u i't / ouubLt 7“ x y r rrssTTTT rc rrs^T . 

JDOF ( 9Qq, 6 ) , £UL (9^,3) , 


data k g p . k 3 * kj,k^/ 


XYZ19DD»3) , 

I WQKQS < ?49 i" 


TrsTSDl - ! jrttt 50 ) 


' SC. 3,900, 

DATA EPS/1. E-IC/ 


M/ 


DATA IvTaWT 
6HFLUI D 


DTR/,Ol7*i532925/» Iz/Q/t 13/3/, 
13 NAME* p IBLANK 


data n j t , notTs* 7*7 


.'JAMES , NAMESt i 
6HGRAVTY,6H5URFTN,6HK1 


, 6H 


I 6/6 

IRTN / 
t 6 H RETURN/ 


1 3 
I H 


TV 
l 6 


IT 
1 8 


IT - 

20 


2 1 
22 


7T 

2*4 


TV 

26 


TT* 

30 


3T 
32 
TT 
3 M 


1 DO 2 
1 30b 


FORMAT 


F o RTTaT - 


t i ox , m 
I 1 OX , 2E 1 7 #3) 


j c i o format 


\ a 6 r 

( 5X , H I 5 ) 


r;zi rwmrr rvrrzTTrn 

102 2 FORMAT (3(5X.EI?«3M 


202 2 


~F ORHAT ( 9 IS) 

FORMAT ( 3 ( 5 X , ! P E 1 >?• . 4 ) 1 


DATA GENERATOR FOR l Oft 60^0 CONTRACT 

CONUINEr. RTTTTU 


xxT s T^ r MTTFTTnr rcuTu ctrrrnrTTiTTi frprn; jtstxt 

9Z DEGREE MODEL- USER SUPPLIED GRID. 

gen e ra t es n ) s tzes and jbof , x~yz , eul matrices on mum for argument 

input to SUBROUTINE FINEL, 

~ZT~ ” ~T " “ ( 2 ) MASS TYPE7 S“TIf _ ”TYPE, DENSITY, ' BU'Ck MODULUS , GRAVITY, 

28 C FLUID ELEMENT JOINT NUMBERS ON N U T E L TO BE READ IN 


t suflKour ines Fluid, graVTy; 

C. SYMMETRIC, A i\ T I ““SYMMETRIC CASE, THAT IS, 

T U«[)X«bUMtI HT HG , V^oYsSOmEThInG", »»DZ*D ON XT PLANE, 

C U*l)X*5i V = DZ-C, wa-OYeSOMETHlNG ON XZ PLANE. 

T X= L 0 N 'MINER AXIS OF SYMMETRY T+UP T , Ys+RjGHT, Z»+INTQ PAPER o 
C EULER ANGLES ONLY USED ON BOUNDARIES WHERE CONSTRAINTS ARE APPLIED , 



54 C READ FLUID SURFACE GRID POINT NUMBERS FrqM AXIS OUT. 


S, M.MFS, I , K G r ) 










N I G « ID P 0 ! N T NUMBERS AT STATE MEN 


CULA ! E XTT. JDW7T.UI MATRICES. 
VAX I s = X TT2 { I ,2 ) 

SECT a iiSIECT ” 

T X I H C = 9 :./5ECT 


0 2 5 I » 1 , K J 

2 5 J D 0 F U » J ) = '0 • * 

CALTTTC' aTET N u M a e s 0 f JOINTS On X - AXJS* 
npax * Z 

“~TTd ”'5 ?■ I 6 P -T,"N G P — 

If I Aast XYT 2 ( 1 GP» 2 )-YAX f S) .LT« EPS) nPAXbNPAXM 


L 

X I S JOIN TT>. T Z "5 at JOINT I » TZ* 3 AT OTHER JOINTS 
QO l "2 I 6 P s 1 t N P A X 
XYZIlHP.n » XVTZm»i r ,T) 

X Y Z ( Iti P I 2 ) = X Y T Z ( 1 G P , 2 3 


I'n i j y ( n 1 i 
NpAXPl = hPAX+1 


l 2H Do 129 IPl.AN£s*,NPLA N £ 


IX = r L U A T ( I rLANE - l ) * T A [ i\J C 
C T X = CUS(TX*DT«| 

“ st x * STn ( r x W ' 

YLOCAU - XYTZ ac,P,2) -YAXIS 


J 3 J+ ] 

XYZU, I) = XTT2(iGP»l) 


loop = I D OF-*- i 


J o a r t j , i ) * loop 

IF (iFCrt . EQ • 1) GO TO 1 2 7 

I D Cfr 5 I 0 OF -*■ I 
J 0 0 F l J » 2 ) = 1DOF 


I OOF * 


Jl> OF ( J » 3 J » IDOp 
129 CONTINUE 



NUTXYZ) (IcUL ( I » J ) , I » 1 » N J ) , J= 1 « 3 I 


REAL JOINT NUMBERS OF GRID POINTS* 

IGP “ 3 


U O Z05 J*1 , N J 

Ip (A&S<XYZ(J,3)~ZAXIS) *GT. EPS ) GO TO 235 


IGP = I G r + I 
JRN (IGP) s J 


*35 LON I JNUt 

CALL iVRlTlM (JRN* NGP , 1 3HJRN • KG P > 


C 

C calculate element JOINT numbers for subroutine FLUID 


write i nutel > nzu na m em,name< » iolank » i&lank 



f?) IGP1 , IGP2, IGP3* I GPR 

If iigpj , e a * a ) so ro 269 


. 3} X E L s 9 


A 

S ON X-AXIS* 
1 a iTTPi ” 

JR * IGP2 



” t- u - i* c. u » 

22s write t nutel *221 q ) nel » j t . J2 » j3 8 jm , 1 z , 1 1 » 1 z » \ 1 


1 














160 

05 

* J ft N U G P 3 ) + I 

se ct- 1 


yxx~ 

03 

* 0 2 + I 



170 

0 6 

* J 5 + 1 



1 7 1 

‘ HEX 

= N £ L + | 



172 

r'7-3 

235 WRITE ( i 1 U T E L » 2 0 l 

•0 1 MEL , J! .02 ,03 ,J4 , J5 ,J6 , J Z . t Z 


ELEMENT is hot on x-axis, 

5- TF' [ & t, L •£<?, HT"GO TO 
DO 2 5S iSECTsl , ivlSFCT 

Ul - JK.j(j“GFT ) "M SE.CT-J 

02 s JrtN < I GP2 ) + 1 SECT- l 


i 1 u r J ) T 1 3 1 1 , l * ! 

JM = U I + J 

3g - , J 2 + l 

J6 = J3tt 

~ NFL = " hFl +1 

25S WRITE l i J U T E L i 2 T I > ) N d L i J 1 . J 2 , J 3 , J *4 , J 5 , J 6 , I Z , ! Z 


= JKUUGP3I + ISECT-1 
a JR/M IIGP^U+I SEC r-.t 


05 = Ol+l 

06 - J 2 + 1 
0 7 = J 3 + 1 
08 * JH +) 

NEL * "N eTm “ 

2&S WRITE. (NOTE L»2-1J£> N E L > J t , 02 . J 2 , J H , 0 5 , 0 6 , J 7 . J 8 


269 WRITE 

tWUTEL .2 

510 

1 1 . r z , 1 z , iz * iz , 1 1 , 1 z , 1 z , 

I Z 

C4LCULA f E 

ELEMENT 

JOINT NUMBERS F 0 R SUBROUTINE 

G ft A V T r 

Aft 1 T a 
•V R I 7 u 

j] .vrrr - 

< HU TEL » 1 
1 NU.TEL . 1 

as 5 ) 

0 2 1) 

NAMES 

I3LANK,NAMEK 



I <NUTEL » 2*22 I GX,GY»$Z 

A 

N 


H < IGP2)+ISECT-t 
= 02+1 

- .MEL+1 ” ‘ 

27S write (mUTEL»2wJ5 ) N E I. , J I . J 2 , J 3 » \ 2 
DO 2 85" I 5 = 3 , ft F S 
I GP » = IFSUS-li 





El * NEL+I 
2 B 5 WRITE ( N U T E L r » 2 C 1 C I N t L , J \ , 0 2 , J 3 , J 4 


write <NUTEL,2C15I 1 1 * 1 1 , 1 2 , I 2 , I 2 




«SSS5* 







M W N 












PH1LIPBIN207*FI.LBFIME 


r 

2 

T 

M 

k 

T 

8 

9 

10 
1 1 
1 2 
1 3 
1 H 

Ts“ 

1 6 
T7“ 

1 e 

T9 

20 

2 1 
22 
TT~ 

2*4 
2TT' 
26 
“27“ 
2 8 

30 

TT 

32 

TT 

3 *♦ 

TIT 

36 

T7~ 

38 

T9” 

HO 

Trr 

'4 2 


CUMP 1 LER ( XM* n rTtWIV^CTHTl " 

subroutine lbfine inutel.nutxyZ,nutm,nutk,nutlt , nutst ,nutmx,nutkx^ 

* ' rvTUT 3 x , N u T 1 * i^TTT 2 t N ut 3 ) 

c 

T MAIN HKU^Ra.TTC!' ( X V 2 ) ", ( JDOF ) , (EUL) AND CALCULATE {ON OPflu^i) 

C ASSEMBLED FINITE ELEMENT MASS, STIFFNESS MATRICES, 

CALLS fWft' WBROUT I NES fIN£L8»YIN , YvsrITE. 

C DEVELOPED BY ft BENF1EL0, C BODLEY, R PHILIPPUS, R WOHLEN. JULY 1973a 
£ LAST REvTSTlTT&'V R”l ftOHLEN* F E 0 R TAR Y 1 9 7 5 . 

C - 

DOUBLE PRECISION V 

COMMON / DOUBLE / V ( 1 2010 ) , LV<12C50) 

" DIMENSION X Y2 ( 2?C? ■ 3 ) , JOOF ( 2CC?,6 ) » £UL(2CC0«3) 

EQUIVALENCE 1 X Y l ( \ ) , V t l) ) , 1 E Vl (1 ) > V ( ) ) 1 ( JDOF(l) ,LV ( 1 ) > 

dTTa ”kr~x t nrcTf7“ "KRJ , KCJ, KRE, KCE, KV / 

* 2-:o:» 3, 6 , 2000, 3 , ! 2 300 / 

T R'TAD STZ'jUOO'r.EiTL nijTXYZ CREATED In BaIA GENERATOR 1 8 D 5 E N • ~ 

REMIND- NUTXYZ 

READ { N U T X Y Z I N J , N C X Tn R 3 , NC J * N« E , NCT 

NERR0k=1 

■' " IF (NCX . NE. 3) GO '"TO" ?9 9 ~ ” — 

NERR0RK2 

If ( NR J 7W7 FTT ,dR. KJ . NE. GO To 99? 

NERR0R*3 

IT TTIHfc. . NE • NJ .ow,” MCE » NE * 3) GO TO W 9 ~ “ ” " ” ‘ 

NERR0R*4 

“ ip i nj.gt.k** Vor. ncx .gt.koTTqr . ~ “ ~ ” ~ 

* NRJ.GT.KRJ .OR. NCJ.GT.KCJ .OP. 

* N h L * 0 I « < R E . 0 R ♦ iT CTTSTTirnrT 5 0 To ~777 

read (nutxyzi ( ( jdof ( 1 , j) , i*i ,mrji , J«i ,ncji 
~ ” read ( n u T x y 2 > 'T ( xyz“([,ui > i*i,nj) , j * f, nc x ) 

READ (NUTXYZ) (( EUL< I i J) , 1*1 » NRE > » J = 1 ,NCEI 

Call fin elb ( x y z , j d 0 f , e u l”, n ut el.nj, 

* NUTM , NUTK , 

* v,lv,kv,krx,<rj,kTe, _ 

* WUTMX , NUTKX , nut 1 , NUT 2 , NUT 3 1 

■ Call Y w r i t E (NUt m. mhmaSs rvTTV ,kv ) ” 

call yvuri te (nutk,mhstif , v ,lv ,kv ) 

— return * " “ " ' 

c 

— v v 9 caul zTwrm rzrruwrim , 
end 



(5PRT FJ.MOTITL 






PHIL1PBIN207«F1**IDTITL 


T 

2 


roMPirf^^rxwiT^ (equi v*cmn7~ — — ~ 

subroutine mdtitl < pt i tle .mode ,freq , nrwt ,kpti tl ) 


c 

C SUBROUTINE to FORM MODE NUMBER and FREQUENCY TITLES for PLQT3. 
Z D t V b. ETTP L u bT^a - s'en'f UXoT 


FEBRUARY 1974. ~ 

LAST REVISION BY R A PMlLlPPUS* MARCH 1975* 


DIMENSION PT I TLE ( KPT ITL > 

z 


10 

1301 

FORMAT ( 3 A 1 0 ) 





1 1 

1 00 2' 

FORMAT C & A 6 > 





1 2 

2 0 3 1 

FORMAT (HHMODE,lM 

»6H, 

F ss»F10*6»6H 

HZ* ) 


I 3 

C 






I'M 


REMIND N R sii T 





t 5 


WRITE {nR'AT , 2051 ) 

MODE 

.FREQ 



I 6 


REWIND NRWT 




• 

1 7 


IF (KPTlfL.EQ. S“) 

READ 

I NRWT f 1 001 ) 

( PTl TLE ( I > 

, l * 1 . 3 ) 

1 8 


IF <KPTiTL.EQ*l 3) 

READ 

( NRWT , 1022 > 

( PTl TLE ( I ) 

,1*1,5) 

1 9 


RETURN 





23 


End 







I3PRT FI* MOD ED 



PHItIP9lM207*FUMQOED 


1 COMPILER ( X M = 1 ) , I E Q U rv*CW! 

2 SUBROUTINE M 0 D E D INUTM l NUTK l NUTP t «UTp*NRSVT! l NUTaitNUTIxn 

T C” “ “ ' ' 

4 c SUBROUTINE TO COMPUTE MOOES USING sparse mass aNO stif matrices 


? C W l ! H afiTSE MODE SUBROUTINE, 

6 C DEVELOPED q Y Vi/ a atNF I eld . M a Y 1978. 


r — 

a 

C“ 

c 

L AS 1 REVISION 

S Y R L ' WO-H L E N . FEBRUARY 1 9 7 S • 


9 “ 

10 


double precision v 

common / dousle / Miis.its), sui5,u5), 2 < 1 1 5 > . «<ns). 


n 

12 

c 

1 


i l 1 5 ) » V< 3300) , LV(35C3) • INORDSI 205) 


l 3 
1 8 

c 

OaT A 

K A , K V 

7 US, 3-3 3 3 / 


IS 
1 6 


Call 

call 

~Y 5 T 0 0 

ystoo 

<nutW|A*nra,nca*ka * k a , v v l v , k v , n u t a i > 

(NUT* , S.NRMS ,NCS ,<A , K A , V , L V , K V , N U T 3 1 ) 


1 7 

i a 


call 

call 

MODE 1 
■VR I TE 

<A,S,-V2*vV,FRE3,NRMS l '>,3,KA jNUTRI) 

. ( &2.NRM5 , .1 %zmz 1 KA ) 


i 9 
20 


call 

call 

WRITE 

WRITE 

(FREQ»NRMS,1 ,8HFREQ*KA) 
IA.NRMS, NR MS.SH MODES, KA) 


2 1 
22 

c 

CONVERT 

call 

DENSE 

YDTOS 

TO SPARSE FOR PLOTS* 

(FREQiNUTF iNRMS* 1 ,KA , t ,v ,LV .KV.NUTSl ) 


23 

28 


call YDTOS 
If ( nrsvt i 

IA f NUTP # PiRMS»NRHS#KA|KA*V i LV,KV,NUT8n 
• LE • 0) RETURN 


ZS 

26 


REWIND NRSVT 1 

call atape t #2 ,nrm5i i ,2h#2 »ka .nrsvtu 


27 

28 


call 

call 

AT A PE. 
Ya T A P E 

tFRE^»NRMS» 1 »<IHFREQ*KA , N R S V T 1) 
( NUTP , SH MOOES i V,LV,KV t NRSVTI ) 


29 

30 


call itaPe 
return 

t N R S V T 1 > 


31 


EnO 





rap r t 


F 1 .rtODES 




20 C 

2 } C OLMNi r iUN U F 1 in HUT VARIABLES* — — 

22 C N = (NUMBER OF MOOES WANTED. 

71 C Wtf = NUMBER Of K A Y L E t G H - R I T Z MOOES TO U5E, 

2q C SHIFT = SHIFT VALUE To USE* 

75 — — — c — maxi i « haxtwu 'm Number of 'iterations to "be performed* 

2 6 C 

27 DOUBLE precision V 

28 COMMON / DOUBLE / V(1$92C)» LV<15 , ’2'3), 7£>» *1 t 7 3 ) » F R E Q < 7S»i 

TV 9 I LORDS' H? 3? ) 

30 dimension m 70, 70>, st 7?, 7C» 

31 C 

32 EQUIVALENCE < VI36AI ),S( H), t L V ( 3 6 4 I ) , A ( 1 ) > 

“33 C 

3*t OaTANIT,nOT/5»6 7 















56 call yarite t nutz , sHflooes » v ,lv,kv j 


57 

58 

-t 

-CONVERT OtlMSt TO bPnKSE I- OR PUOTT^ - 

call ydtcs t freo # nutf ,nu 1 1 , ka » i , v ,i v , kv ,nuT7 i 


59 

60 


If (NRSv-Tl .L£. di return"' 

REWIND NRSVT1 


51 


LaLL" ft! A PE ( 4 2 »NU , t , 2 M IV 2 ,KA ,NK5VT U 


62 


call wtape ifreq ,nu , i . hhfreq ,<a ,nrsvt i ) 




call y*tape inutz ,5hmooes , v ,i_v ,kv ,nssvt d 


64 


call ltape < nrsvt \ ) 


65 

C 



66 


return 
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Tno 




IsFFT rTTSTTTST 





ITE ELEMENT . » . 

: STIFFNESS MATRIX 


r F OK A SURFACE TENSION TRIANGLE ELEMENT -ViTri UNRE5TRA ined au 
C STIFFNESS MATRIX IS IN (SLO 3 A L COORDINATE DIRECTIONS* 



W AKE TRANSLATIONS, 

EULEr ANGLE CONVENTION is GL03AL X.V»Z P£ RftfUT A T I ON* 
CALL'S hURiT'A" SUBROUTINES fT A 3 A . DC OS 2 , K 2S T 1 , l ZB 0M3^ 
DEVELOPED 3 Y Rl *0rtL£N« FEBRUARY 1 975c 



SUBROUTINE ARGUMENTS 


» INPUT MATRIX OF GLOBAL X , Y , Z COORDINATES AT TRIANGLE 
ROWS 1 » 2 * 3 CORRESPOND TO X,Y»Z COORDINATES* 

COLS 1,2*3 CORRESPOND TO JO I NTT - ! , 2 » 3 • SIZE (3 ,3) . 

* INPUT MATRIX OF EULER ANGLES (DEGREES) AT TRIANGLE JOINTS 
— ROaS ),2“,"3 CORRESPOND TO GLOBAL xTyTZ PERMUTATION* 

COLS 1,2,3 CORRESPOND Tq JOINTS 1,2,3* SIZE(3,3). 


= input surface tension, cforce/length) , 

E k = INPUT TYPE OF STIF M A T R I X W A N T g D « _ 


, USES K2ST 1 , LINEAR DISPLACEMENT FIEL 
a OUTPUT STIFFNESS MATRIX * S I Z E ( <? , 9 ) . 


PACE MATRIX. SIZE! 18,1°) • 
input row dimension of cj in calling program* 




Input row dimension of s in calling program* min*?* 
rNPTJT ROW DIMENSION OF - 3Ti IN CALLING PROGRAM* MIN=~i 



If i k s 


. L r * 9 .or 


SORT ( ( C J ( 1 


K W \ • 


C J ( 1 , l M **2 + (CJ(2,2)|?CJ(2, l ) >»*2 

+ ( C J ( 3 , 2 >“"C J ( 3 , n ) «»2 ) 


CJ< l ,2) )**2 + (CJ(2,3)"CJ(2,2>>*«2 

+ (CJ(3|3)-CJI3,2)»*»2) 

S NRTT~IC J (|,3)-CJ(1,IM*»2 + (C J ( 2 , 3 ) -C J( 2,1)) **2 

+ (Cjf 3t3>'~CJ(3» 1 > W«2> 


* ( S L 1 3 * * 2 + 5 L 1 2**2-SL23**2) / ( 2*3*SLI2 1 

* SNRTISL 13**2-X3**2 J 

( N A MEK. *EN- 6HK1 ) GO TO lit! 



2 5 T l a LINEAR DISPLACEMENT FIELD* 
no call K2su isli2,X3, y3,st,s,ks) 


call DC0S2 ( C J , E J , W I , K C J « K E J . K«V t } 


SELECT 01 ROWS, 


UQ Z 10 J= 1 , 9 
DO 213 1*1 ,3 


( 1 1 j 1 * 

DO 2 15 J*1 ,3 















i«t c STIFFNESS HAI'UX is in global coordinate DIRECTIONS* 



2 4 


ROtfS 1,2,3 CORRESPOND To X,Y,Z COORDINATES, 




26 

' M 

mm 

input 

COLS I , 2 , 3, 4- CORRESPOND TO JOIN » 5 I, 2 ,3,4. SUtUill. 
MATRIX OF EULER ANGLES (DEGREES ) aT QUAO JOINTS. 

27 

c~ 



"WttS 1,2,3 CORRESPOND TO GLOBAL X,Y,Z PERMUTATION* 

28 

c 



COlS 1,2, 3, 4 CORRESPOND TO JOINTS I, 2, 3, 4. SIZEt3,4). 

TV 



sn 

"Tnfut 

SURFACE TENSION (T ORCE /LENGTH) 

30 

c 

N A rs E K - 

Input 

TYPE OF ST IF MATRIX v! ANTED. 

3 1 

c~ 



- KI, USES 4 TRIANGLES, OVERLAP AVERAGE, 

3 2 

c 

S = 

output 

STIFFNESS MATRIX 5 I Z E C (2,12). 


33 (T W7 “ INPUT 

3 4 C ,«»2 » INPUT 


vV ij KK SPACE MATRIX. SIZE 1-9 , 9l *' 
WORKSPACE MATRIX- SIZE! 14, 1«) 




I NPUT RO,< DIMENSION OF CJ IN CALLING PROGRAM* MIN»3* 
Input row dimension of ej in calling program* min«3« 
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5 6 

2 30 

E ,V II , 3 ) 

S EJI1,3) 



S7 


call si f 2 $t (cr*";F?r,sT ,namek i »* 2 * 3 , 3 , k a 1 ,k *2 ) 



SB 


call rev add ( .s,*i » i v \ , 1 v 1 »s» 9*7*12,12* kwi,ks) 



59 


Do 221 I 

= 1,3 



63 


C *V ( I , 1 ) 

S CJU.l) 



6 I 


t # M . 1 ) 

= EJ( I ,1 ) 



62 


C y< (1*2) 

* CJ(1,3) 



63 


Eft C X v 2 1 

= EJ ( I ,3) 



69 


C W ( I » 3 ) 

= C J ( I , 4 ) 



6 S 

201 

E iv { 1 , 3 ) 

= F. J ( I , 4 ) 



6 6 


call stf 

2ST ( CV» , E* ,ST ,M ArtEK ,W 1 t 2 , 3 , 3 ,K1S 1 ,KW2 ) 



67 


call rev add < . $ , * 1 » 1 v , 1 v 2 > s , 9,9,12*12 , k»i,ksi 



68 


Do 203 I 

*1,3 



69 


~Ot (l»n 

- C J ( 1 , 1 ) 



70 


E -H ( I » 1 1 

= Eju,n 



7 1 


C # ( I , 2 ) 

« Cjl I ,2) 



72 


E ft C I * 2 ) 

= E J U ,2) 



73 


t W (1,3) 

- C J ( 1 ,4 > 



79 

203 

EW( I *3 ) 

= l J (I , 4 ) 



. 75 


call s t f 2 5 t * c * , e vd , s t , n a * e k ♦ w 1 » w 2 , 3 , 3 , k « 1 , k w2 j 



76 


call rev ado (•5«<ni*iv3,iv3*s. 9,9,12,1.2, kwi,kS) 



77 


DO 2-35 1 

-1,3 



7 0 


c an i , n 

= C J ( I ,2) 



v 7 9 


t w < r, i ) 

= E J ( I ,2 ) 



80 


cw( 1 ,2) 

= C J ( I ,3 > 



8 1 


Eft f I ,2 ) 

= EU ( I ,3 ) ~~ “ 



82 


C W ( I , 3 ) 

= C J ( I , 4 ) 



83 

235 

EW ( 1 ”, 3 ) 

s E J ( I ,4 > 



89 

\,x 

Call STK2ST ( C %■ , £ W , ST , N A UEK , iV l , ,S2 , 3 * 3 , Ktf 1 , Kfl2 > 



BIT 


CALL - RLVAO'O < . 5 , W 1 , ( V 9 , IVT, S* 9 , 9 , ! 2 , i 2 , K # 1 , K S ) 



86 


RETURN 




87 

* c 





88 

999 

call zz&omb ( ahstpsst , nerror > 



89 


end 
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P H I L 1 P B l '■J 2 U 7 * F" l • S U R F T N _ 

j. — — iL£rt"7 xm« i ) ,‘i e .jutwcthti 

2 SUBROUTINE. SURFTN ( xrz , JDOF ,EUL jNUTEl *nj, 

^ — ■■ - T TTU TKX j v , T , S » K X , K J , K E , K JV ) 

n 0 i M £ N S I 0 fo X r 2 ( K X ♦ ! ) * J D 0 F ( i< J , I ) » E U L ( K E t l ) , ( K VV f n , T ( K jy t 1 ) I S t K IV » t ) 

5- u * h t msi a n cju»4),r:j[3^),ivin2) 

6 DATA NIT.NOT/ S,6 / 

— ______ A T A 'NXMEL/6H5URTTn7, 13LWK/6H " /. KCJ/3/" 
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2M 
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C 
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C 

~C 

C 

T* 
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T 
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“C 
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5 U 3 W 0 J f I N E TO" CALCULATE' ( ON OPTION) FINITE ELEMENT” , e , 

STIFFNESS MATRICES AND IVECS (ON N I J T K X ) t 

FOR S 1 J R F A C E TENSION ELEMENTS* TRIANGULA* (JOINT M = 0) OR 

QUAD!?! LATERAL (JOINT M »GT. 1). ______ 

stiffness matrices are in global coordinate DIRECTIONS# 

global COORDINATE OR DER is 

( u » V , W J JOINT T7~ THEN JOINT 2 , 3 , { M > . 

H HERE U.V,* ARE TRANSLATIONS® 

I VEC GIVES ELEMENT DOF INTO -GLOBAL DOF, EXAMPLES... ~ 

IVECC6)s83M PLACES ELEMENT DOF 6 INTO GLOBAL DOF 3 3_M c 

"T VTC ( 3 j ONITS ELEMENT OOF 3 FROM GLOBAL OOFo THIS CONSTRAINS 

i ELEMENT i> 0 F 3 TO Z ERO MOTION® 

DAlA ARRANGEMTNT ON N'J TKX FOR EACH FINITE ELEMENT IS <W = K) 

Uh-UTE (NUT A' X ) NAHE/V*NCL t N R,NC*NAMEL i I IBLNK , I s 1 »5> , 

((rtU,J),| = l,fJR),J=!,«C)|(IVEC(I),I = l,NC) ” 
CALLS FORMa SUBROUTINES P AG E HD , S TF 2 5 T » S T F 3 S T , 2 ZB 0 M S , 

DEVELOPED 3 Y R\ -^OhLEN. FEBRUARY 19757 


ft * ft * * * <e # a V o ft* ft * * * ft » * * * ♦ ft * 


« ft ft 


INPUT DATA REAO IN THIS SUBROUTINE FROM NUTEL* IF NUTEL * NIT* OATa I 
W LA U I- HUM CARUSO — ““ 


Na-MEK FORMAT (A6) 

S T : format isx»eid> 

' 20 NFL, J\ I J2* J3, JM FORMAT < 5 I 5 y 

IF (J1 ,E0* 01 RETURN 

GO TO 20 


DEFINITION. OF INPUT VARIABLES. 

NAMF.K w rv^e - OF STI^f-NGSS MATRIX WANTED.' 

= K I , LINEAR DISPLACEMENT ASSUMED* 

«— 03 6 HN 0 STIF, NIT - stiffness mat^Tx calculFIXIt; 

ST * SURFACE TENSION (FORCE/LENGTH)* 

TTFC = n NITTTLeMENT n umber* FWTFFTrENCE 0NlT7 NO I USED IN 

CALCULATIONS. WRITTEN ON NUT<X* 

tti s joint number ai Element vertex i. 

J 2 = JOINT NUMBER AT. ELEMENT VERTEX * f 

TT" = JOINT NUMBER At - ELEMENT VERTEX 3 * 

jm 5 joint number at element vertex m. (used for quadrilateral). 

'I HE" Cl E n L N I M a V BE NUMbE RcO cTTOoT^T^ tTf? C 0 OTT'TTR - CLOCKWISE. 


EXFLANATT O- N U T — PTPUT FORMATS® NUMB'VR indicates Card COLUMNS UsEO^ 

l * INTEGER data, right ADJUSTED. 


DECIMAL FTmn data, ANYWhERT-TTj" FIELD 
CARD COLUMNS SKIPPED* 


(I4t*,***i*,.»*ft*« *"• ft***aftft**o«*** ft ft * «~~ 


E X P 0 N E N IT” R I G hi I AD JCTSTETi? 


55 


r 


SUBROUTINE ARGUMENT S“TaT:l'“T HP U TT 













MATRIX of joint GLOBAL X»Y ,Z LOCATIONS. HOWS CORRESPOND 
TU~J~CHTN T fTCTHB E K S . C 0 L U H N S ~T,2 .3 CUR RESPOND To f HE J u I nT 
x,y,z locations respectively. sizeujj.3). 

RaTT? I X 0 F JO I NT GLOBAL UEGREES iJF“FREEDOM. kO W 5 COKKtSP 
TO JOINT N 0 M i. j E i< S • COLUMNS 1.2.3 CORRESPOND TO THE JOINT 


C NOVEL 


CORRESPOND 
i - THE JOINT 


TO THE JOINT 



6 C 

ROTATION OOFS. S I t E ( H J ♦ 6 > • 

H A T R IX "D F JOINT E U L ETT N'GLES (DEGREES). ROWS CORRESPOND 
TO JOINT NUMBERS* COLUMNS 1,2*3 CORRESPOND TO THE 
GLO'BAL X , Y * Z PERMUTATION^' SIZE(mJ,3J* 

LOGICAL NUMBER OF Tape containing element input data for 


SUBROUTINE. IE N U T E L * N I T , DATA IS READ FROM 
NUMBER of JOINTS OR ROWS IN MATRICES <XYZ)» (JDOF) 
LOGIC AL NAJM&E.R OF UT 1 Cl T Y _ Ta'PE ()N WHICH ELEMENT 
STIFFNESS MATRICES aNP ivecs ARE OUTPUT. 

UTrnCX'T‘.<uT T;E ZERO I F STTTFNESS MATRIX IS NOT FORME 
USES FORTRAN READ. NR HE- 


MATRIX WORK SPACE, m l N SIZE! 12, \2), 

MATRIX WORK SPACE. M I N S I l E ( 9, 9} . 

MATRIX WORK "SPACE . MIN S I 2E ( I 8 , i 8 ) « 

ROW DIMENSION OF XYZ IN CALLING PROGRAM. ' 

ROW 0 I ME NS ION OF JDOF IN CALL I N G PROGRAM* 

ROW DIMENSION OF' EUL IN CALLING PROGRAM. 


ROW DIMENSION OF N , T, AND S IN CALLING PROGRAM. HIN=18. 


ICC I FORMAT” 
IZZ 2 FORMAT 
I 333 FORMAT" 
2 SD 1 FORMAT 


2532 Format <//2?x 


2C.?3 FORMAT 


203 R format 


( A 6 1 

( 5 X t E 1 G ! . ^ ) ' 

tsisF 

1//2SX HIM INPUT DATA FOR SURFACE TENSION STIFFNESS 


6H< triangle orquadr i lateral ) elements) 
i//2?x 8jn input data for surface tension stiffness 

■ Hsu { TR I ANGLE or QUADRILATERAL) ELEMENTS (CONTIWUI 
1/12X7HSTIF ■ A 6 , / 1 &X8HST * IPE10.8, 

// l 5 X 7 H EL E ME NT 13X7 HJOl NT 1 tiX7HJ0INT 2 J 3 X 7 H J 0 I NT 3 
13X7HJ0INTH 


/ 1 5X6HNUM3ER > 
(18X»5(l5fISX)) 


N LINE = 


AUEHD 
WR I TE (NOT ,202 1 ) 


READ (NUT£L,13Q1) N A M E n 

READ (NUTEL, 1 D 0 2 ) ST 

Hfrt I TE ( N 0 T , 2 0 D 3 ) MAMEKtST 

READ (NUTEL , I0Q3 ) N E L , J 1~, j 2 ,jT7JM 
IF (J1 .LE. 55) RETURN 


) NAMLK.ST 


G T , N J .OR 


Nj * 0 R • J3*GT#NJ *0R* JH#GT#NJ) G 






na 



1 12 

c 



TT 3 — * 

~c 

h U N I't ^ i w n t LLtnt^l CUUK[) I A T fc, LULA! 1 u N b , u v l_ b 'X ANOLLb, KcVAUU IVtC. 

1 M 



DO M2 I *1 ,3 

' IIS 



T3 rrrn = m ztj rrn 

1 l A 



c j 1 1 » 2 ) = XYZuz.n 

TT7 



LJlltJi = X T £ ! J 3 » 1 1 

\ \ s 



EJ( 1 il 1 = EUL ( Jl , I ) 

1 l 9 



“TjT 1,7) - E’jLIJ 2, 1.1 

J zo 



EJ ( I ,3 ) « EUL ( J 3 , 1 ) 

n?i 



1 V I n ) " = 'JO Of ( J l , I ) 

122 



1 V 1 ( 1*3) = JOOF (vl2, 1 1 

\ 2 3 



1 VI N+o) = J D 0 f { J 3 » I ) 

I 2 4 



IF <04 .C.T. Q) GO TO 4 4 

rzb 




l 26 




177 



350 ? 0 u 

123 

r 



rzy 


■ 4*f 

“DTI *T5 T=T 71 ' ” 

1 3Q 



C J ( 1 » 4 } = X Y l ( 0 4 , n 

rrr 



E- kj l i 1 ) * t U 1 J S 1 I 1 

132 


4 S 

1 V 1 l l + 9 ) = JDOF { J 4 , I ) 

” 133 



n c 0 l - rz 

13 m 



call stf3st ( c J ,e j ,st , w , t , s ,kc j ,kc j ,kw »k .« ,kv< > 

TTS 

— r 



1 36 


l J c 

If 1 fc A M E K • E Q * 6 H .OR. N A M E K .E°t 6 H N 0 S T IF ) GO TO 20 

TT 7 



' iTrRTT0"R"s7 

138 



IF {NUTKX * L L • t?) GO TO 999 - 

” TTV 



W RI f L TmJTOI NAHEK * >'! E L , N COL , NCO L » M A M F L » < I B L W K » I = \ ,5 > » 

I 43 



* t C Ml ( 1 ,J 1 »I«l ,MC 0 L> iJ-t ,NCOL) t ( 1 V l I 1 ) » 1 * 1 » ncol ) 

nn 

C 



1 ^2 



GO TO 27 

i *T3 

77 



144 


999 

call zzbohb i ahsurftn .nessor » 




End 


wir 
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Input requirements for the Free Surface Static Equilibrium Shape 
Computer Program are given in Section 6.2.1. Input requirements for 
the Vibration Analysis Computer Program are given in Section 6.2.2. 


6*2.1 Input Requirements - Free Surface Static Equilibrium Sha pe 
Program - A description of the input requirements to the free surface 
static equilibrium shape program along with listings of sample input 
are presented in this section. The first three cards required are to 
satisfy subroutine START (as explained in Reference (3)). Data input 
to either the search mode or the survey mode is accomplished by means 
of the NAMELIST facility available in Univac 1108 Fortran V. Section 
6.4 of Reference (6) contains a detailed explanation of this facility. 
Input values which must be specified for the survey mode are: 

AC0F0 the base value for the A sweep, i.e., the first A value 

used will be ACOFO+DACOF 


BONDNO the nondimens ional ratio of inertial forces to surface 

tension forces, based on container length 

DACOF the increment applied to AC0F0 to generate successive values 

of A 


DELTAT the arc length increment used in the numerical, integra tion 

algorithm to generate a solution trajectory 

NA the number of A values to be computed, i.e., NA trajectories 

will be generated ranging from ACOFO+DACOF to AC0F0+NA*DAC0F 

IPRNT trajectory results from the numerical integration will be 

printed every IPRNT integration intervals, ignored if PRINT= 
.FALSE. 


NX the number of tank axis intercepts used to generate solutions 

from a given trajectory (A value), i.e., for a given value of 
A NX solutions will be generated ranging from XUP-DX to XLO+DX 
where DX= (XUP- XLO) / NX 

PRINT controls the printing of intermediate results, PRINT=.TRUE. 

prints all intermediate results, PRINT=. FALSE . prints only 
the solution summary (Figure 3-1) 

RMAX the maximum radius of the container 

SEARCH controls the selection of program mode, SEARCH^. FALSE, for 

the survey mode 
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XLO specifies the lower limit of the tank axis intercept for which 

solutions are generated from a given trajectory 

XMAX length of the container measured along the axis of symmetry 

Xup specifies the upper limit of the tank axis intercept for 

which solutions are generated from a given trajectory 

Figure 6-2 shows a listing of a sample problem input data for the survey 
mode. Note that BONDNO is not given in the input data, thus the default 
value specified in the computer program (Section 6.1.1) is used. 

Input values which must be specified for the search mode are: 


ACOFO 

BONDNO 

DACOF 


DELTAI 


initial value of A 

same description as survey mode 

initial increment to be applied to ACOFO, the program tries 
both ACOFO+DACOF and ACOFO- DACOF in searching for an improved 
solution; if none is found, DACOF is halved and the search 
repeated. The value of A corresponding to the desired solu- 
tion must lie in the range ACOFCH;DACOF for this procedure to 
work 

same description as survey mode 


EPSC 


value of the error function at which convergence is esta- 
blished 


IPRNT 


same description as survey mode 


NX same description as survey mode 

PHID the desired value of contact angle in degrees 


PRINT 

RMAX 

SEARCH 

ULPCT 

XLO 

XMAX 

XUP 


same description as survey mode 

same description as survey mode 

same description as survey mode, however, SEARCH= . TRUE . 

for the search mode 

desired value of ullage volume percentage 

same description as survey mode 

same description as survey mode 

same description as survey mode 
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Figure 6-3 shows a listing of a sample problem input data for the search 
mode. Note that SEARCH is not given in the input data, thus the default 
value specified in the computer program (Section 6.1.1) is used. 

Multiple runs maybe made in either mode by repeating the cards 
required by subroutine START and the NAMELIST data as many times as de- 
sired. The run is terminated when START reads the word STOP in the run 
number field of the first card. 
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SVYTUG 

SAMPLE SUHVEY MODE RUN FOR 
TUG-LIKE TANK. 

$ T NO AT A 

ACOEO-O . 1 ) oooo • 
t) ACOF = - 0 ♦ 0 1 0D0O . 
DELTAT=0. 20000* 

NA= 1 5 • 

MX=40 . 

PRTMT= .FALSE. > 
RMAX=32 . 2000 • 
search^ .false . . 

XLO-O.OD0O. 

*MAX=1 46.5DO0 • 
xup=] 46,^000 * 

SEND 

STOP 





Figure 6-2 SAMPLE PROBLEM INPUT DATA - FREE SURFACE STATIC 
EQUILIBRIUM SHAPE PROGRAM, SURVEY MODE 



BONDTUG WARNS* 

TUG-LTKE CONTAINER f BOND NUMBER = l.o* CONTACT ANGLt = 0.3 
ULLAGE VOLUME = 80.0 PCT 
IINDATA 

ACOF0=0.061562BDC0, 

BONO wo = 1 - odoO . 

DACOF = 0.000031 P5D0 0 * 

DELTAT = 0.£OT>00« 

EPSC-0 * BO-05 * 

NX = 50 * 

PHT 0=0 .3000 ♦ 

PRINT*. TRUE . * 

PMA X = 3.2 • ?D0 0 i 
ULPCT=8.0. 0O0 0* 

XLG=A6 . ODOO « 

XMaX=T4'6.5D0"0V ‘ 

X U P = A 8 . 0 * 

SEND 

STOP 




Figure 6-3 SAMPLE PROBLEM INPUT DATA - FREE SURFACE STATIC 
EQUILIBRIUM SHAPE PROGRAM, SEARCH MODE 
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6.2.2 In p ut Requirements - Vibration Analysis Program - An explana- 
tion of the input to the vibration analysis program, along with a listing 
of input data to a sample problem (see Figure 6-4) are given here. Input 
formats to subroutines START, READ and READIM are explained in Reference 
(3). 


Card No. 

(Ref. 

Figure 6-4) Input 


1 

Run no . , 

cols. 1-6; name, cols. 11-28 

2 

Title 1, 

cols. 1-78 

3 

Title 2, 

cols. 1-78 

4 

1 INI TIL * 

or 'NOINIT 1 , cols. 1-6 


Three cards to satisfy subrou- 
tine"START" 

To initialize or not to initialize 
the reserve tape 


5 

6 

7-30 

31 

32 

33 

34 

35 

36 

37 

38 

39 

40 


Number of sectors in 90° model 

Matrix name (XYTZ), no. of grid 
points, no. of cols. (3). 

Grid point X, Y coordinates and 
0 Z values 
Ten zeros 

Z coordinate of longitudinal axis 


Format (10X,I5) 

l 


Subroutine READ 


Format (10X,E10) 


Grid point of corner if container 
is a cylinder 

Matrix name (GP-CW) , no. of rows 
(1) , no. of grid points on con- 
tainer wall 

Container wall grid point numbers 
Ten zeros 

Matrix name (GP-FS) , no . of rows 
(1), no. of grid points on fluid 
surface 

Fluid surface grid point numbers 
JJen zeros 

Mass option; Ml for lumped, M2 
for consistent 


Format (10X,I5) 


Subroutine READIM 


Subroutine READIM 


Format (A2) 


41 

Fluid mass density, scale factor x 
bulk modulus, surface tension 

Format 

(3(5X,E10)) 

42 

Acceleration 

Format 

(5X,E10) 

43-57 

Element number, grid point numbers 
for elements (CW numbering) 

Format 

(415) 

58 

Ten zeros 



59 

Mode calculation option; MODES for 
large sparse, MODED for small dense 

Format 

(A5) 
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60 

No . of modes wanted 

Format 

(10X,I5) 

61 

No. of modes used 

Format 

(10X,I5) 

62 

Shift value for (convergence 

will be about this value) 

Format 

(10X,E17:0) 

63 

No. of maximum iteration allowed 
(cards 60-63 are omitted if MODED 
option used) 

Format 

(10X, 15) 

64 

Plot option; PLOT or NOPLOT (if 
NOPLOT , omit cards 65-87) 

Format 

(A6) 

65 

First mode plotted (usually 1) 

Format 

(10X,I5) 

66 

Last mode plotted (<no. of modes 
used) 

Format 

(10X,I5) 

67 

No. of plot views 

Forma t 

(10X, 15) 

68 

Stereo plot option (1 for stereo, 
0 if not) 

Format 

(10X,I5) 

69 

No. of tracings of line 

Format 

(10X,I5) 

70 

Print option (1 - yes, 0 = no) 

Format 

(10X, 15) 

71 

Plot title 

Format 

(13A6) 

72 

Elements plotted; GRAVTY plots 
surface joints, FLUID plots all 
joints 

Format 

(A6) 

73 

No. of view positions 

Format 

(10X,I5) 

74 

Roll angle 

Format 

(10X,E10) 

75 

Matrix name (C0EL0C) , no. of rows 
(1) , no. of cols . (3) 



76 

X, Y, Z coordinates of center of 

\ Subroutine READ 


eyes 



77 

Ten zeros 



78 

Matrix name (VPLOC) , no. of rows 



79 

(1), no. of columns (3) 

X, Y, Z coordinates of viewpoint 

r 

\ Subroutine READ 

t 

80 

Ten zeros 



81 

Read plot data option (1 - yes, 
0 = no) 

Format 

(10X,I5) 

82 

Cross-section option (1 = yes, 
0 - no) 

Forma t 

(10X,I5) 

83 

Read modal data option (1 - yes, 
0 = no) 

Format 

(1QX,I5) 

84 

iS ^^tual mode number of first mode 
^Mcula ted 

Format 

(10X,I5) 

85 

Scale factor on modal displacements 

Format 

(10X,E10) 
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86 

Option on superposition of undefonned 
and deformed joints (l^es, 0=no) 

Format 

(10X,I5) 

87 

Symmetry option. Use XSA in cols. 
15-17 



88 

Plot option; PLOT or NOPLOT (if 
PLOT repeat cards 65-87) 

Format 

(A6) 

89 

End of data . STOP in cols. 1-4 




/ 
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Figure 6-4 SAMPLE PROBLEM INPUT DATA - VIBRATION ANALYSIS PROGRAM 
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($ IF 

(A NT| 


41 NU 

g SHIFT 
MAXIT 
44 PLOT 
4S MN START 
(A MN END 

41 NVIFWS 
STEP 
NTRACE 
70 JF PRIMT 
11 IUS OX TANK 
GRAVTY 

73 N VIEW RT 
74 roll ancle 

15 COFLOC 1 

7fe 1 1 

•n oooooooooo 
78 VPLOC 1 

71 i 1 

oooonoonoo 
| IF READ PD 
CROSS SECT 
IF read MD 
REAL mnI 

scalf fact 

IF SU p R p OS 
SYM 

NOPLOT 
STOP 


6 

3,0E-6 
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Figure 6-4 SAMPLE PROBLEM INPUT DATA VIBRATION ANALYSIS PROGRAM (cont'd) 



